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ABSTRACT 

Recent interest in pollution control and the proximity of Monterey 
Bay to the Naval Postgraduate School prompted an investigation of the 
circulation in the bay* The first phase of the study consists of 
solving the simple cavity flow problem. A vorticity-streara function 
relationship is solved using an explicit, time dependent, finite dif¬ 
ference sclierae. Solutions for selected Reynolds^ numbers and length 
to width ratios of the cavity are obtained. Values are chosen to give 
an indication of the flow patterns occurring over a wide range of 
these parameters. 

Equations for a refined model are derived to include the effects of 
the bottom topography, frictional forces and the Coriolis force. A 
numerical procedure similar to the one applied to the simple cavity 
flow problem is used on the refined equations. The topography of 
Monterey Bay is used in this study. Results indicate that closed 
circulation patterns are more probable in Monterey Bay than in other 
emba>Tiients due to the presence of the submarine canyon. 


2 








TABLE OF CONTENTS 


I. INTRODUCTION.^5 

A. BACKGROUND.’ . . . . ^5 

B. NUMERICAL MODELING.^9 

C. SCOPE OF PRESENT PAPER.20 

II. NUMERICAL METHODS .2.2 

A. INTRODUCTION.22 

B. NUMERICAL RELATIONSHIPS ..'-2 

1. Taylor’s Theorem . 22 

2. Representation of Derivatd.ves .23 

3. Representation of Arakawa Jacobian . 25 

C. RELAXATION METHODS.27 

1. Residuals .27 

2. Relaxation.28 

3. Overrelaxation.29 

4. Optimum Overrelaxation Factor, Ropt . 29 

a. Introduction.29 

b. Theoretically Determined Value . 29 

c. Empirically Determined Value . 

D. CONSERVATIVE AND TRANSPORTATIVE PROPERTIES . 34 

E. CONVERGENCE AND STABILITY .38 

III. GOVERNING EQUATIONS . 40 

A. BASIC EQUATIONS OF MOTION . 40 

1. Horizontal Momentum Equations . 40 

2. Continuity Equation . 41 

3. Integrated Equations of Motion . 42 


3 





















































B. DEVELOPMENT OF VORTICITY EQUATION . 43 

1. Derivation .43 

2. Nondimensionalization . 45 

3. Nondimensional Parameters . 50 

C. DEVELOPMENT OF STREAM FUNCTION EQUATION . 52 

1. Derivation .52 

2. Nondimensionalization . 53 

IV. SIMPLE CAVITY FLOW .54 

A. GOVERNING EQUATIONS . 54 

B. NUMERICAL APPROACH .. 56 

1. Upwind Differencing for Solving the Vorticity 

Equation .56 

2. Method for Solving Vorticity Equation Employing 

the Arakawa Jacobian .57 

3. Optimum Overrelaxation to Solve Poisson's 

Equation for the Stream Function . 58 

4. Boundary Conditions . 59 

a. Inflow.60 

(1) Velocity Profile . bO 

(2) Stream Function Profile . b3 

(3) Vorticity.66 

b. Walls.67 

c. Downstream Continuation . 68 

d. Lid .68 

e. Corners.69 

5. Convergence and Stability Criteria . 70 

6. Computational Sequence . 71 

7. Programming Optimization . 72 

C. RESULTS.73 


4 































D. CONCLUSIONS 


78 


V. REFINED MODEL . 79 

A. GOVERNING EQUATIONS . 79 

1. Bottom Friction Equations . 79 

Final Form for the Vorticity and Stream Function ... 81 

3. Convergence and Stability.U5 

B. INVESTIGATION OF THE RELATIVE IMPORTANCE OF THE 

CORIOLIS FORCE AND BOTTOM FRICTION . 86 

1. Coriolis Force . 86 

2. Bottom Friction.89 

C. IRJMERICAL APPROACH .90 

1. Method for Solving the Vorticity Equation . 90 

2. Optimum Overrelaxation for Solving the 

Stream Function Equation . 90 

3. Boundary Conditions . 91 

a. Inflow.91 

b. Lateral Boundaries . 93 

c. Downstream Continuation . 95 

d. Lid .95 

e. Bottom Topography . 96 

4. Computational Sequence . 98 

D. RESULTS. 98 


1. Neglecting the Coriolis Force and Bottom Friction 


2. Neglecting Bottom Friction . 

E. CONCLUSIONS.115 

APPENDIX A - Development of the Bottom Friction.HI 

APPENDIX B - Solutions for the Simple Cavity Flow Problem . 123 

REFERENCES .1-^8 


5 





























INITIAL DISTRIBUTION LIST . 1A9 

FORM DD 1473 .153 


6 






% 

I 


f 

I 

i 


5 


\ 




I 

I. 


■ 





LIST OF FIGURJilS 


Figure Page 

1.1 Location of Monterey Bay . 17 

1.2 Monterey Bay . 18 

2.1 Numerical grid . 31 

2.2 Plot of overrelaxation factor vs. number of 
passes required for convergence for Re = 1, 

10, 50, and 100 and aspect ratio of 1 . 35 

2.3 Plot of overrelaxation factor vs. number of 
passes required for convergence for Re = 1, 

10, 50, and 100 and aspect ratio of 2 . 36 

2.4 Plot of overrelaxation factor vs, number of 
passes required for convergence for Re = 1, 

10, 50, and 100 and aspect ratio of 3 . 37 

4.1 Specification of boundaries. 60 

4.2 Velocity boundary layer profiles . 61 

4.3 Unsuitable boundary layer profile . 64 

4.4 Stream function boundary layer profiles . 65 

4.5 Portion of numerical grid plotted. 74 

4.6 Gomparison of four solutions for the stream 

function. 76 

4.7 Sample vorticity plot, Re = 10, B = 2. 77 

5.1 Frictional coefficient, R1 . 83 

5.2 Frictional coefficient, R2 . 84 

5.3 Grid for refined model . 94 

5.4 Hypothetical bottom topography . 97 

5.5 Bottom topography simulating Monterey Bay . 99 

5.6 Bottom topography with no smoothing passes . 100 

5.7 Bottom topography with one smoothing pass.101 


7 


























I 



m 




5.8 Bottom topography with two smoothing passes . 102 

5.9 Bottom topography with five smoothing passes .... 103 

5.10 Refined model, hypothetical bottom topography, 

volume transport stream function, Re = .01-100, cF = 0 105 

5.11 Refined model, topography simulating Monterey 
Bay, volume transport stream function. 

Re = .01-100 cF = 0.106 

5.12 Simple cavity flow, stream function, Re = .01, B = 2 , 108 

5.13 Refined model, hypothetical bottom topography, 
volume transport stream function, Re = .01-1.0> 

eF = .23 no 

5.14 Refined model, hypothetical bottom topography, 
volume transport stream function, Re = .01-1.0? 

cF = 23.0.Ill 

5.15 Refined model, bottom topography simulating 
Monterey Bay, volume transport stream function> 

Re = .01, eF = .23 x 10“^ .112 

5.16 Refined model, bottom topography simulating 
Monterey Bay, volume transport stream function, 

Re = .01, eF = .23 X 10^ .113 

5.17 Refined model, bottom topography simulating 
Monterey Bay, volume transport stream function. 

Re = .01-1.0, eF = 1.0 x 10^.114 

B.l Simple cavity flow, stream fimction. 

Re = 1, B = 1.124 

B.2 Simple cavity flow, vorticity, Re = 1, B = 1 .... 125 

B.3 Simple cavity flow, stream function, 

Re = l,B = 2...126 

B.4 Simple cavity flow, vorticity, Re = 1, B = 2 .... 127 

B.5 Simple cavity flow, stream function, 

Re=l,B=3 . 128 

B.6 Simple cavity flow, vorticity, Re = 1, B = 3 .... 129 

B.7 Simple cavity flow, stream function. 

Re = 10, B = 1 .130 

B.8 Simple cavity flow, vorticity. Re = 10, B = 1 .... 131 


8 













I 

I 

I 

I 

{■ 

I 

I 





B.9 Simple cavity flow, stream function, 

Re = 10, B = 2 .132 

B.IO Simple cavity flow, vorticity, Re = 10, B = 2 .... 133 

B.ll Simple cavity flow, stream function. 

Re = 10, B = 3 .. 134 

B.12 Simple cavity flow, vorticity. Re = 10, B = 3 .... 133 

B.13 Simple cavity flow, stream function. 

Re = 100, B = 1 .136 

B.14 Simple cavity flow, vorticity. Re = 100, B = 1 ... 137 

B.15 Simple cavity flow, stream function. 

Re = 100, B = 2 .,.133 

B.16 Simple cavity flow, vorticity, Re = 100, B = 2 ... 139 

B.17 Simple cavity flow, stream function. 

Re = 100, B = 3 .1^'5 

B.18 Simple cavity flow, vorticity. Re = 100, B = 3 ... 1^1 

B.19 Simple cavity flow, stream function, 

Re = 1000, B = 1 .l'^2 

B.20 Simple cavity flow, vorticity. Re = 1000, B = 1 ... l'^3 

B.21 Simple cavity flow, stream function. Re = 1000, . . . 144 

B = 2 

B.22 Simple cavity flow, vorticity. Re = 1000, B = 2 ... 1^3 

B.23 Simple cavity flow, stream function. 

Re = 1000, B = 3 .l^'^S 

B.24 Simple cavity flow, vorticity. Re = 1000, B = 3 ... 1^^ 


9 












I 




LIST OF TABLES 


Table 

I 

II 


Page 

Tabulation of the rate of convergence of Ropt .... 33 

Frictional coefficients . ^2 


10 










f 









<r ^ 


TABLE OF SYMBOLS AND ABBREVIATIONS 


a 


B 

D 

e 

f 

F 

g 

h 

H 

i 

J 

In 

log 

L 

0 

P 

r 

Rl, R2 

Re 

Ropt 

t 

T 


Argument in Ekman solution for slope currents 
Coefficient of horizontal eddy diffusivity 
Coefficient of vertical eddy diffusivity 
Aspect ratio, AX/AY 
Characteristic depth 
2.7182818- 

Coriolis parameter, 2 ^ sin 
Nondimensional coriolis parameter 
Acceleration due to gravity 
Depth 

Nondimensional depth 

Jacobian 

Natural logarithm 
Common logarithm 
Characteristic length 
Order of magnitude of 
Pressure 
Radius of earth 

Frictional coefficients 

U L 

00 

Reynolds number, 

Optimum overrelaxation factor 
Time 

Nondimensional time 


11 



i 




Tr 

Transport 

u,v 

x,y components of velocity 

u,v 

x,y components of nondimensional velocity 

U 

00 

Characteristic velocity 

W 

u + iv 

x,y ,z 

Space coordinates 

X,Y 

Nondimensional space coordinates 

z 

Nondimensional vorticity 

a 

/2/2 

3 

Beta plane approximation to coriolis force 

t 

Time increment 

t 

c 

Critical time increment 

Ax, Ay 

Grid spacing 

V 

Del operator 

V 

Nondimensional del operator 


Laplacian 


Nondimensional Laplacian 

3 

Partial differentiation 

e 

BL^ 

Nondimensional parameter, ^ 

oo 

C 

Vorticity 

K 

Von Karman’s constant 

T1 

Height of sea surface 

9 

o 

Latitute 

TT 

3.1415926- 

P 

Density 

bx by 

T , T 

x,y components of bottom stress 

WX \7Y 

T , T 

x,y components of wind stress 


12 











'i' 

Subscripts 

x,y 

w 

Superscript: 

k 

n 


Turbulent shearing stress 
Stream function 

Nondimensional stream function 
Angular velocity of the earth 


x>y grid positions 

Partial differentiation with respect to x,y 
Wall value 


Iteration level 
Time step level 


13 


I 

I 


i 


I 

f 






ACKNO^^^EDGEMENTS 


Dr. Edv/ard B. Thornton proposed the topic for this thesis and 
guided me during my investigation. I am especially indebted to him 
for his assistance, encouragement and professional advice. 

Several other members of the Naval Postgraduate School staff were 
helpful in the development of this study. Dr, Jerry A, Galt assisted 
me in the physical interpretation of the problem and also suggested 
the approach used to incorporate the bottom friction into the refined 
model. Dr, Henry Crew and Dr. Frank D. Faulkner aided me in the pro¬ 
gramming and mathematical portions of this thesis. 

This study would have been impossible without the facilities and 
cooperation of the staff at the William R. Church Computer Center, 
which is part of the Naval Postgraduate School, All programming was 
done for the IBM 360/67 Computer, which is located at this computer 
center. 


14 






I 



I 


1 



I. INTRODUCTION 


A, BACKGROUND 

Monteiey Bay, located on the West Coast of the United States 
approximately one hundred and twenty miles south of San Francisco, 
presents a good opportunity to study currents within an embayment 
(Figure 1.1). This bay is vjell suited for circulation studies because 
of the following geographic and oceanographic features. The bay is a 
large one and, to a close approximation, it has symmetrical boundaries. 
It is semi-elliptical in shape with a major axis of twenty miles and 
a width of about eight miles (Figure 1.2). Bisecting the bay is a 
large submarine canyon. This canyon originates in the bay offshore of 
Moss Landing. At the mouth of the bay, eight miles from Moss Landing, 
the canyon attains a depth of six-thousand feet. This topographic 
feature is as large as the Grand Canyon in vertical relief. The effect 
of the Monterey Submarine Canyon on the circulation within the bay is 
one of the problems which is investigated in this study. A desirable 
oceanographic feature is the presence of currents which flow past the 
bay. The Davidson Current flows north along the West Coast of North 
America from November to February. Flowing southward during the 
spring and summer is the California Current. It is assumed in this 
study that oceanic currents, such as these, are the driving force for 
the circulation within an embayment. The currents drive the circulation 
by advecting momentum into the area in question. The effects of tidal 
forces and local winds are neglected. It is not meant to be implied 
that these forces are unimportant, but that only a single component of 
the total system is being considered in this analysis. 
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Monterey Bay is conveniently located for conducting circulation 
studies for several reasons. First, the Naval Postgraduate School is 
located on the shores of the bay in the city of Monterey. Both the 
Postgraduate School and the United States Navy are interested in circu¬ 
lation studies of this kind. The facilities at this institution are 
well adapted to carry out such studies. The Department of Oceanography 
at the Naval Postgraduate School is interested in circulation studies 
and has access to oceanographic research vessels. Also, located at the 
School is a large computer center with a well-trained staff and an 
IBM 360/67 Computer along with all the necessary peripheral equipment. 
Secondly, there are many communities located on the shores of Monterey 
Bay which are becoming increasingly aware of pollution problems. An 
important question that needs to be answered before pollution can be 
controlled is: what is the effect of currents on the discharged 
effluent? To answer this question the circulation regime of the area 
involved must be determined. It is hoped that the results of this study 
will provide some answers to the above question. 

Prior to this study, little had been done to ascertain the circu¬ 
lation patterns within Monterey Bay. This thesis is part of a larger 
research effort being conducted by the Naval Postgraduate School aimed 
at determining the forces which dominate in controlling current patterns 
within an embayment and the types of patterns which are generated. Pre¬ 
sently, there are three phases to the project. Two phases concern 
numerical modeling. One numerical model, the subject of this thesis, 
assumes oceanic currents drive the circulation within a bay. A second 
model simulates the tide and wind induced circulation. The third phase 
involves field studies conducted in order to assess the models which 
have been developed. 
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B, NUMERICAL MODELING 


Recently, with new developments in partial differential equation 
theory and finite difference methods, the field of fluid mechanics has 
progressed tremendously. The advent of large-scale computers has 
spurred these developments and made it possible to solve partial dif¬ 
ferential equations numerically which have proven, so far, to be 
impossible to solve analytically. Due to the great speed of these 
computers it has become possible to solve large systems of simultaneous 
equations rapidly enough to put numerical prediction models in fields 
such as Meteorology and Oceanography on a real time basis. 

Most circulation problems, basically, involve the solution of a 
large number of simultaneous equations \7ith specified boundary con¬ 
ditions. Practically all studies of this type are concerned with 
either small-scale flows with simple geometries such as: channel flow, 
flow around a cylinder, and simple cavity flow or with large-scale oceanic 
circulation studies. Little has been done to try to apply these equations 
and numerical techniques to medium-scale flow such as the circulation 
within an embayment. There are several possible reasons why more medium- 
scale flow problems have not been investigated. One reason is the high 
cost of studies of this type. Computer time is very expensive and not 
always readily available. Medium-scale problems are in general highly 
variable and each case must be dealt with individually. Another reason 
is that, before the increased awareness of pollution problems, there was 
little interest in the regional circulation beyond that caused by waves 
and tides. Lastly, the application of the equations of motion to large 
or medium-scale flow problems requires assumptions about the forces and 
terms involved which have not yet been universally accepted as being valid. 
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There are many numerical procedures that can be used to try to solve 
a system of simultaneous partial differential equations. An interesting 
result of applying different numerical techniques to the same problem is 
the possibility of obtaining different solutions. One of the problems of 
numerical Tiiodeling is determining whether the correct solution has been 
reached, if only one exists. There has been a lot of research into the 
problem of the solution of partial differential equations by Ames [Ref. 2]. 
Considering the Navier-Stokes Equations, he proposes that a unique 
solution exists only below a certain Reynolds’ number and that above this 
Reynolds’ number several solutions exist; whereas, above a larger critical 
number no solutions exist. The critical Reynolds’ number marks the 
transition to turbulent flow. To illustrate the existence of multiple 
solutions to a partial differential equation, he gives an example of ci 
quasi-linear, elliptical partial differential equation for which there 
is no unique solution. 

An extensive study of various numerical procedures used in the 
solution of compressible and incompressible laminar separated flow for 
simple geometries was made by Roach [Ref. 7]. His dissertation is quite 
thorough in its consideration of the problems encountered and the 
techniques used to resolve them. Some of Roach’s techniques are 
extended in this study to the problem of motion within an embayment in 
the simple cavity flow section of this thesis. 

C. SCOPE OF PRESENT PAPER 

This thesis attempts to take the equations of motion and apply them 
to the problem of the circulation in an embayment and solve them 
numerically. A progression is made from the simple cavity flow problem 
to the refined model. The simple cavity flow problem considers: local 
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rate of change of vorticity with time, advective terms, and lateral 
shear stresses. In addition, the refined model considers: planetary 
vorticity tendency, topographic vorticity tendency and bottom friction. 
Procedures for numerically solving these equations are examined while 
economy of computer time and the validity of the solution are kept in 
mind. An attempt is made to generalize the computer program so that it 
can be used on any embayment. The generality is obtained by allowing the 
input of an arbitrary bottom topography. 

The model is tested on Monterey Bay in order to discover, not only 
the circulation patterns resulting from the interaction of the forces 
involved, but also the relative importance of the forces. To determine 
the importance of the forces involved, various computer runs are made 
with different combinations of the forces acting, and the effect on tlie 
circulation patterns is examined. It is hoped that the results of these 
investigat;Lons will shed some light on the processes which control circu¬ 
lation patterns and help provide answers to some of the problems which 
now exist. 
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II. NUMERICAL METHODS 


A. INTRODUCTION 

This section consists of a brief discussion of some of the aspects 
of numerical modeling which are used in this thesis. There are many 
excellent texts which contain thorough discussions of this material; 
therefore, the material is not covered in depth. All of the repre¬ 
sentations derived in the following section, II-B, will be used in the 
development of the numerical equations solved in later sections. .The 
following paragraphs introduce these concepts and develop the necessary 
relationships. 

B. NUMERICAL RELATIONSHIPS 

To represent a partial differential equation in numerical form the 
starting point is usually Taylor’s Theorem. From this theorem, the 
necessary expressions which represent derivatives in their various forms 
can be derived. Numerical relationships for terms such as the Laplacian 
and the Jacobian can then be developed from these expressions. Basically, 
these expressions determine the value of a derivative at a given point 
in terms of the values at surrounding points and the grid spacing. 

1. Taylor’s Theorem 

The finite difference representation for the derivatives of a 
function at a given point in terms of the values of the function at 
surrounding points can be derived from Taylor’s Theorem. 

Taylor’s Theorem in Two Dimensional Space: 

If u and its derivatives of order ^ p+1 are single valued, 
finite and continuous at every point on the line segment 
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from (x,y) to (x+Ax,y+Ay) then there is a point 
(a,b) on that line segment such that 

P 

u(x+Ax,y+Ay) = u(x,y) + I ~~ N ^ ^^.l) 

(N+M)=l 9x 9 y^ 


where the remainder, R, is given by 


,N+M 


R = 


^ u(a,b) 

^ •' M (^y) 


N+M^p+l 9x^9y^ 


It can be seen that if Ax and Ay < 1 as N and M -> «> then R ^ 0; 
and in the limit, the expansion equals the function. If y is held 
constant (i.e., Ay = 0) there is a contribution to the summation terms 
only when M = 0, therefore 


u(x+Ax,y) 


u(x,y) 


P 

+ E 
N=1 


n: 


d U 


Bx 


/A ^ / P+lx 

(Ax) + 0 (x ) 


( 2 . 2 ) 


If X is held constant, (Ax = 0), similarly 

P M 

u(x,y+Ay) = u(x,y) + I ^(Ay)^ + 0 (y^"^^) (2.3) 

M=1 ‘ By 

Using these equations, expressions for derivatives at a point 
are developed in Section II-B-2. 

2* Representation of Derivatives 

The following grid illustrates, diagrammatically, the 
notation that is used to fix the grid position of a given point. 


23 










I 




I 


I 


x-Ax,y+Ay 


x,y+Ay 


x+Ax,y+Ay 


^ ^ 


/ 



""i-lj+l 

-J.! 

1,1+1 

Vl,j+1 

Ay 


x-Ax,y 

\ 

x,y 

x,y+Ay 

e ^ 

• 

0 


"i,j 

""1+1,1 

x-Ax,y-Ay 

x,y-Ay 

x-hAx,y-Ay 

• 

• 

• 

u. . . - 

1-1,1-1 

“i,j-l 

""i+i,i-i 


Using Taylor’s Theorem with y held constant Equation 2,2 
the new notation becomes 

3u. . 


u. 


.2 

. 3 u. . ^ ^ 

•J .1 • = Ax + —-(Ax) + 0 (Ax) 

1 + 1 ,j i,j 8x 2 


If - Ax is substituted for Ax the equation becomes 

8u 


u. 1 . = u. 


, a^u. , 

1 _llJ. 


. , . Ax + -^-(Ax)^ + 0 (Ax)^ 

1 - 1,1 1,1 3 x 2 3^2 


Linear combinations of 2.4 and 2.5 with the higher order terms 
neglected yields 

Representation Mag, of Error Diff. Scheme 


u. . - . - u. . 


3x 

J J__ 

Ax 

9u 

U. , - u . 

1,1 1-1,1 

9x 

Ax 

9u _ 

u .,1 . - u. , , 

1+1,1 1-1,1 


9x 

.2 


2Ax 


0 (Ax) 

0 (Ax) 

0 (Ax)' 


Forward 


Backward 


Central 


u. . + u. - . 

3 u ^ iH-l,i 1-1 ».1 


- 2u. 


9x 


(Ax)^ 


0 (Ax)' 


in 


(2.4) 


(2.5) 
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Similarly, derivatives with respect to Y can be found. 



Representation 

Mag, of Error 

Diff, Scheme 

3u 

u. - U. , 

1,3+1 1,3. 

0 (Ay) 

Forward 

9y 

Ay 

3u 

U. . - u. 

1,J-1 

0 (Ay) 

Backward 

h 

Ay 

9u 

'^i.i+l ■ '^i,3-l 

0 (Ay)2 

Central 

9y 

2Ay 

9^0 

u. +u. .,-2u. . 

1,3+1 1,3-1 

0 (Ay)^ 


a 2 
9y 

(Ay)^ 

. 


A lepresentation for 

2 

the Laplacian, V u, 

can be developed 


from the above expressions, 

2 2 

2 _ d u d u 

V u = - TT + —;r 


ry u.,- . + u, - . - 2u. . u. . + u. . - - 2u. . 

^ 1+1,J 1-1,J i,J. ^ 1>J+1 1,J-1 1,J 


( 2 . 6 ) 


(Ax)‘ 


(Ay)' 


3, ReT)resentation of the Arakawa Jacobian 

This study utilizes Arakawa’s technique for representing the 
Jacobian which is an average of three different expressions for the 
Jacobian [Ref, 3], The method Arakawa developed is used for two 
reasons. First, it has the desirable property of conserving certain 
quantities two of which are the vorticity and the vorticity squared, 

A further discussion of the conservative property is contained in 
Section II-D, Secondly, the representation is a very stable one which 
helps to minimize convergence problems in solving the numerical equations. 
The price that is paid for this stability is the high degree of smoothing 
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inherent in this representation which may eliminate some small-scale 
features. Although these features may be real, it is not felt that 
their elimination is very important because this study is only inter¬ 
ested in the general flow patterns. Examination of small-scale features 
in the flow patterns with the simplifying assumptions which are made in 
this study would not be realistic. 

In order to develop the desired relationship for the Arakawa 
Jacobian, the following three expressions for the Jacobian are 
considered, 


J(A,B) 


9A 3B _ 3A 8B ^ A 
9x 9y 9y 9x 


(2.7) 


J(A,B) = 


J(A,B) 


^ (A^) 

9x ^ dy^ 


9y 9x^ 


^ (A-^) 
9y '■ 9x'' 


9x 9y^ 


= J 


= J 


( 2 . 8 ) 

(2.9) 


The Arakawa Jacobian is given by 


. i jB ^ 


Using the same notation as in the previous section, the expressions 
become 


t(A.^n .-A. , .)(B, .^,-B. . ,) - (A. . ,)(B.,t .-B. ) 

4AxAy 1+1,J 1-1, J i,J+l i,J-l i,J+l i»J-l 1+1,J 1-1»J 


J® =tA- tA,^, ,(B,^, ,_,) - A,_, , (B,_, ,_.) 


4AxAy i+l,j i+l,j+l i+l,j-l i-l,j i-l,j+l i-l,j-l 


- A, .j.i(B.,t .,1-B, , .,,) + A, , ,(B.,^ . ^-B. . . .)] 

i,J+l 1+1,j+1 1-1,j+1 i,J-l 1+1,J“1 1-1,J-1 


[B, (A, 


- B, , , (A, , ,^,-A. , , ,) 


4AxAy i,j+l i+l,j+l i+l,j-l i,j-l i-l,j+l i-l,j-l 


-B.,t ,(A.,- ,,--A,,, . ,)+B. , .(A. T ..^-A. T • i)] 

1+1,j 1+1,j+1 1+1,j-1 1-1,J 1-1,J+1 1-1,J-1 
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Collecting terms and factoring yields the desired form. 


J(A,B) = 


. (A, , ,+A,^, , ,-A, + 


12AxAy i+l,j i,j-l i+l,j-l i,j+l i+l,j+l 


B . (A_ -^A. ^ •j.i“A. ^ ,-A. ^ ..i) + 
i,j+l 1+1,J 1+1,j+1 1-1,J 1-1,j+1 


B. ^ . (A. ..,+A. ^ -.I “A. . i”A. . -) + 
1-1,j i,j+l 1-1,j+1 i,j-l 1-1,j-1 


B. , T (A, . ,+A, . , 1-A,^. .-A.., . ,) + 
i,j-l 1-1,J 1-1,j-1 1+1,J 1+1,j-1 


B-.i -.1 (A.., .-A, , .) + 

1+1,j+1 1+1,J i,J+l 


( 2 . 10 ) 


B. T ,.T (A, .,,-A, , ,) + 

1-1,J+1 i,j+i 1-1,J 


B. ,, , - (A, .-A, , .) + 

i~l,J“l i“l,3 i,J“l 


Bj.i • -I (A. . i-A.., .)] 

i+l,j-l i,j-l 1+1,J 


C. RELAXATION METHODS 

One general approach to solving a set of simultaneous equations is 
the relaxation method. This method attacks the problem by determining 
the difference between an estimated solution to the equations and the 
actual solution. Each successive solution to the equations differs 
from the ’*true” solution by a determined amount. The aim of this method 
is to manipulate the equations so that the differences from the "true” 
solution approach zero. The manner in which this is achieved is 
explained thoroughly in Ref, 1 and discussed briefly in the following 
three sections, (II-C-1, II-C-2, II-C-3). This discussion begins with 
the definition of a residual, 

1, Residual 

One can write a set of simultaneous equations such that all the 
terras are on one side of the equality sign, for example: 
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X + Y - 3 = 0 


Residuals are defined as quantities which represent the values of the 
equations when values are assigned to the variables, X and Y. The 
solution to the equations is reached when the residuals equal zero. 

X + Y - 3 = R (2.11) 

X 

2X + 3Y + 2 = R (2.12) 

y 

The residuals, R and R , are a measure of the closeness of X and Y to 
^ X y 

the exact solution. The problem of solving simultaneous equations 
then becomes one of finding values for X and Y which make the residua]s 
approach zero. One method for doing this is called relaxation. 

2. Relaxation 

The relaxation method involves taking the largest residual and 
computing a value for its associated variable which makes the residual 
zero. Using the above two equations as an example, the variable X is 
arbitrarily associated with the residual in Equation 2.11. Likewise, the 
variable Y is associated with the residual in Equation 2.12. 

Assuming that R is the larger residual, then to make R =0 

X X 

only the value of X would be adjusted. Now that R^ = 0, is the 

largest residual, and the value of Y in Equation 2.12 would be adjusted 

to make R =0. This process is continued until the values of R and R 
y X y 

are within acceptable limits. 

It is not economical in terms of computer time to determine 
which equation has the largest residual and relax only it. Making many 
sweeps through all of the equations and relaxing them simultaneously is 
a more efficient method. The reason for this is that logical statements 
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take much more computer time to execute than arithmetic statements. 

The process which relaxes all of the equations simultaneously until 
all of the residuals are smaller than a predetermined value is called 
simultaneous relaxation. Modifying this procedure, it is found that 
even more computer time is saved if simultaneous overrelaxation is 
used. 

3. Overrelaxation 

Overrelaxation is a process whereby the residuals are not 
reduced to zero but are overcompensated for in order to speed the rate 
of convergence. This is the usual process used because it minimizes 
the computer time needed to solve the equations. This method is 
approached by writing the numerical equations in such a way that the 
residuals are multiplied by a convergence factor which is called the 
optimum overrelaxation factor, Ropt. A critical part of the problem 
is determining the factor which will yieild the fastest convergence rate. 

4. Optimum Overrelaxation Factor, Ropt 

a. Introduction 

A factor, Ropt, can be determined which maximizes the 
convergence rate of a set of simultaneous equations. The convergence 
factor must be between 1.0 and 2.0 and its choice may be critical, 
because small differences in the factor may make large differences in 
the convergence rate. There are two ways of determining Ropt. One 
way is through theoretical considerations and the other is to determine 
it empirically. These methods are discussed in the following two 
sections. 

b. Theoretical Ropt 

Application of theoretical considerations must be restricted 
to simple partial differential equations and simple boundaries because 
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of the complexity of the problem [Ref. 7]. Frankel applied the concept 
of overrelaxation to Leibmann’s method for solving the following 
equation 

Z = (2.13) 

Leibraann’s method for solving the above equation involves deriving its 
numerical representation and solving for 4^. This method produces 


^k+l ^ ^ [4'^ . + ('i'^ + 4'^'*’^ ,) 

i,J 2(1 +b2) 

- (Ax)^ - 2 (1 + B^) 4'’^ .] 

i>J i.J 


. (2.14) 


The bracketed term represents the residual, and the solution is being 
reached as it approaches zero. 

Frankel determined that for a rectangular grid with 
Dirichlet boundary conditions the optimum overrelaxation factor is 
given by 


^ ^ ,1 - v^~a, 

Ropt = 2 [-] 


a 


(2.15) 


where 


and 


B = 


cos (-) + B cos (-) 

a = [-2-^ 

1 + B 


Ay 


N = Number of rows (Ax) 

M = Number of columns (Ay) 
It is found that for a 41 x 80 grid: 
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B = 1 Ropt = 1.8574 


B = 2 Ropt = 1.8916 

B = 3 Ropt = 1.9064 

c. Empirically Determined Value 

It can be seen from the last section that Ropt, theoretically, 
depends upon the grid size ratio as represented by B and the dimensions 
of the array being relaxed. A study was conducted to experimentally 
determine Ropt for the following grid which is the one used in the 
simple cavity flow problem. 


I - 41 GRID POINTS 



Figure 2.1 Numerical grid 

Leibmann’s equation for solving for the stream function is the one 
relaxed. This is Equation 2.14 with Ropt = 1. 

Evaluating the equation involves initializing the grid and 
relaxing the equations until every value on the grid differs from the 
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previous value by less than 1.0 x 10 . This is approximately the 

limit of accuracy on the IBM 360/67 Computer without resorting to 
double-precision accuracy which would si^piificantly Increase compu¬ 
tational time. To further optimize use of computer time the convergence 
is checked only once every ten passes through the grid. Each pass con¬ 
sists of a sweep through the grid from left to right, top to bottom 
followed by a sweep from right to left, bottom to top. A limit of 
two-hundred passes through the grid is placed on each test. 

Results of this study are summarized in Table I and pre¬ 
sented in graphical form in Figures 2.2, 2.3, and 2.4. Table I gives 
the number of passes through the grid required for every value on the 
grid to converge. The criterion for convergence is that the difference 
between two successive passes is less than 1.0 x 10 Results are 
given for four different Reynolds’ numbers and for aspect ratios of 
1, 2, and 3 represented by Ax = .05, .10, and .15 respectively. In all 
cases ^y is held constant at .05. Figures 2.2, 2.3, and 2.4 represent 
the same data in graphical form. This illustrates the effect of Ropt 
on the convergence of the stream function equation in a more explicit 
manner. 

One conclusion which can immediately be drawn from these 
empirical results is that Ropt is also a function of the Reynolds’ 
number. The most significant result of this study is that overestimating 
the value to be used for Ropt may significantly increase the convergence 
rate; whereas, underestimating Ropt is not nearly as critical. A brief 
comparison between the theoretical and empirical results is as follows: 
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Tabulation of Rate of Convergence of Ropt 
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B 

Theoretical Ropt 


Empirical Ropt 




Re=100 

Re=50 

Re=10 

Re=l 

1 

1.8574 

1.650 

1.725 

1.625 


2 

1.8916 

1.875 

1.875 

1.875 

1.700 

3 

1.9064 

1.850 

1.850 

1.825 

1.800 

From the 

i above summary, it can be seen that the theoretical 

value for 

Ropt is 

usually greater than the 

empirical 

value. 

This situation might 

lead to 

a poor choice for Ropt; 

therefore. 

trying 

to apply 

the theoreti 

cal results from simple problems 

to more c 

omplex one might 

lead to a 


bad choice for Ropt. 

D. CONSERVATIVE AND TRANSPORTATIVE PROPERTIES 

One property which is desirable in a finite difference method is 
the conservative property. A finite difference scheme possesses the 
conservative property if it preserves a certain integral relation of 
the momentum equations [Ref. 7]. This integral relation states that 
over some region the time rate of change of a given quantity must equal 
the net flux of the quantity across the boundary plus the production rate 
of the quantity within the region. Whether a finite difference scheme 
preserves the conservative property depends upon the form of the momentum 
equations used and the numerical scheme employed. The possession of 
this property does not imply that the scheme is more accurate than one 
that does not possesses it unless one of the criteria of the problem 
is the conservation of a given property. Upwind differencing, used in 
the simple cavity flow problem, conserves vorticity. The Arakawa 
Jacobian conserves vorticity, vorticity squared, linear momentum and 
kinetic energy; however, all procedures that employ the Arakawa Jacobian 
do not necessarily conserve these properties. 
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2.2 Plot of overrelaxation factor vs. number of passes required 
for convergence for Re = 1, 10, 50, and 100 and aspect ratio 
of 1. 
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2.3 Plot of overrelaxation factor vs. number of passes required 
for convergence for Re = 1, 10, 50, and 100 and aspect ratio 
of 2. 


36 

















NUMBER OF PASSES REQUIRED FOR CONVERGENCE 



Figure 2.4 Plot of overrelaxation factor vs. number of passes required 
for convergence for Re = 1, 10, 50, and 100 and aspect ratio 
of 3. 
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Another desirable property is the transportative property. The 
finite difference form of a flow equation possesses the transportative*. 
property if the effect of a perturbation on a transport property is 
advected only in the direction of the flow [Ref, 7]. An intuitive 
approach to the problem leads one to agree that a perturbation in the 
flow field should only move in the direction of the flow. The method 
of upwind differencing possesses the transportative property. 

E. CONVERGENCE AND STABILITY 

The mathematical basis for stability and convergence lies in .the 
theory of linear partial differential equations. Linear theory is 
used as a guideline for analyzing non-linear equations such as those 
used in this study, 

A convergent finite difference scheme is one in which all values of 
the finite difference solution approach the values of the exact solution 
as the finite difference size approaches zero. All the terms in a finite 
difference representation can approach the corresponding terms in the 
analytic equation, and yet, it is possible that the entire equation will 
not approach the exact solution; therefore, the convergence criterion 
will be unfulfilled. Convergence, therefore, is concerned with the 
limit of the entire equation and not the individual terms. 

Stability is achieved when the cumulative effect of all round-off 
errors is negligible. This implies that the errors do not increase 
exponentially. In the time dependent approach used in both the simple 
cavity flow problem and the refined model, the stability of the numerical 
method employed to solve the vorticity equation depends upon the size of 
the time increment and the degree of convergence of the stream function 


38 


1 ^^ 




equation, A critical time step can be determined for each problem. 
Exceeding this time step causes the equations to become unstable and 
diverge. It is found that the critical time increment depends upon the 
grid spacing, Ax and Ay, the Reynolds’ number and the depth (in the 
refined model), The most stable solutions for the simple cavity flow 
case occur at Reynolds’ numbers of approximately one hundred and stability 
decreases as the grid size ratio of Ax to Ay increases. For the refined 
model, the greatest stability is found around Reynolds’ numbers of one- 
tenth, Two is the only aspect ratio that is used in the refined model. 
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III. GOVERNING EQUATIONS 


A. BASIC EQUATIONS OF MOTION 

The governing equations in this study are derived from the hori¬ 
zontal momentum equations and the continuity equation. Developments 
of the momentum and continuity equations can be found in many books on 
Physical Oceanography such as: The Principles of Physical Oceanography 
[Ref. 5]. Applying Newton’s second law to a continuous volume of fluid, 
the horizontal momentum equations can be derived. The equation of con¬ 
tinuity can be derived from the principle of the conservation of mass. 
Starting with these equations, the complete development of the governing 
equations used in this thesis is derived belox>7. 

1. Horizontal Momentum Equations 

The Eulerian representation of the momentum equations desired 
in this study is derived from Newton’s second law utilizing the following 
simplifying assumptions: 

a. The curvature of the earth is negligible for the distances 
considered in this study; therefore, cartesian coordinates can be used 
for the coordinate system. 

b. Fluid is homogeneous. 

c. Fluid is incompressible. 

d. Hydrostatic approximation. 

e. Vertical component of Coriolis force is negligible. 

f. Turbulent stresses are proportional to the gradient of 
the mean flow. 

These assumptions imply that the only external forces acting are: the 
Coriolis force, gravity, wind stress, and bottom stress. Assuming that 
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primes denote non-integrated variables, the Horizontal Momentum 
Equations are 


9t 


9t 


+ u' 


+ u' 


9u’ . , 

-- h V 

9u' 

- fv' 

1 

9P' 

3x 

9y 

P 

9x 

+ V 

9v' 

+ fu' 

= - i 

9P' 

3x 

9y 

P 

9y 


+ 


^u' + A 


"9z2 


V 


^v' + A 


9 V 


^ 9z2 


(3.1) 

(3.2) 


where in Equation 3.1 


9u' 

9t 


= Local rate of change of velocity with time 


9u' 

9x 


+ v' 


9u' 

9y 


fv’ 

1 9P' 
p 9x 



= Non-linear field accelerations 
= Coriolis term 
= Pressure gradient term 

= Lateral shear stress 

= Vertical shear stress 


2. Continuity Equation 

The Equation of Continuity is needed as one equation in a 
system of three equations in three unknowns. Physically, this equation 
imposes the condition that no two fluid particles can occupy the same 
space at the same time. Continuing the convention that primes denote 
non-integrated variables and neglecting any changes in the vertical, the 
Equation of Continuity can be expressed mathematically as 


3t 9x 3y 


(3.3) 


For Incompressible, homogeneous flow, the density, p, is constant; 
therefore. Equation 3.3 becomes: 
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(3.4) 



3. Integrated Equation of Motion 

Equations 3.1, 3.2, and 3.4 represent three equations in three 
unknowns: u, v, and p. The solution to these equations is fully deter¬ 

mined when proper boundary and initial conditions are specified. Two 
approaches may be used in order to solve these equations for the desired 
integrated form.^ 

a. The equations can be integra1:ed over the entire water 
column, but this cannot be done without the non-linear terms being 
known explicitly as a function of depth. 

b. The real problem can be replaced by one which assumes a 
homogeneous fluid and the last term in Equations 3.1 and 3.2 can be 
replaced by a body force in which u and v are independent of depth. 

The second method is employed; therefore, the last term in both 
equations is replaced by a body force which results from the wind stress 
and bottom stress, defined as 


wx bx r,2 I 

T T d n 


(3.5) 




(3.6) 


Substituting these terms into the horizontal momentum equations and 
dropping the primes to denote that the variables are now integrated over 
the entire water column, the integrated equations of motion are 



(3.7) 


1 . 


The following discussion is taken from Ref. 4, pp 12 
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(3.8) 


9v av 9v 1 9P 2 


ph ph 


Integrating the equation of continuity over depth, and applying 
Leibnitz Rule of integration yields the desired integrated continuity 
equation. If the sea surface is assumed to be of constant height, 
arbitrarily selected as 0, and if h equals the depth of the ocean, the 
equations take the following form 


o 

o 


// 

3v , 
dz = 
3y 

-h 

-h 



Using Leibnitz Rule where u and v are independent of depth 


o o 

J udz + j udz = 0 

-h -h 


The desired equation is 


^(uh) d (vh) 

3x 9y 


(3.S) 


B. DEVELOPMENT OF VORTICITY EQUATION 
1. Derivation 

Starting with the equations of motion (3.7 and 3.8) and the 
integrated continuity equation (3.9), the desired vorticity equation 
with the wind stress assumed negligible is derived below. The local 
wind stress is assumed negligible as this force is not being considered 
because the oceanic current is the only forcing function being examined 
in this study. Assuming the wind stress is negligible Implies 

wx wy ^ 

Cross differentiating the equations of motion, 

3(3.8) _ 3(3.7) 

3x 3y 
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yields: 


8 .8v 8u x ^ . 8v __ 8u x _ 8u x .8v _ 

9t 3x 9y 9x 9x 9y ^ 9x 9x 9y 9y 9x 9y 


, 9 ,9v 9uv , r I i- , 9f 

+ V -r- (-r-■^) + f—+u — + f — + 

9y 9x dy dx 9x 9y dy 

2 9v 9uv p 9 9 ,n 

^9x 9y^ ~ ^3x ^ph ^ “ 8y ^ph ^ ^ 


The vorticity, is defined as 


(3.10) 


- _ .5ii 

^ 9x 9y 


(3.11) 


Substituting into Equation 3.10 and collecting terms gives 


If fc ^ I7 


(^ + ^) 

'‘9x 9y'^ 


(C + f) 


(3.12) 



(I-) _ ^ (I ) 

ph ’ 9x >h ^ 


Expanding the continuity equation produces 


9uh , 9vh _ ,9u , 9v^ 

^ ^ ^ 97^ 


, ^ 9h . 9h 
h + u V — 

9x 9y 


= 0 


which in turn yields: 


9u 

9x 9y 


1 

h 


. 9h ^ 9h. 

<“ t ; 


(3.13) 


Substituting into Equation 3.12 it is found 



U It: (C + f) + V I77 (c + f) - -r (u 


9x 


9y 


9x 


u. 

■^" 97 ^ 


(C + f) = 





(I—) _ ^ ( 1 —) 

Vh '' 9x >h ^ 


(3.14) 


now. 
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“ 3? 37 '“ir> - h “ 3i ‘^ ■'• + h'' -57 « 


1 , 9h I 9h. , , ,, 

;;2 <“ 37 ■"'' 37 > 


(3.15) 


therefore: 

9^ . , r 9 ,^+f, . 9 ,^+fs, . _2 .9 ,T^^s 9 

9t-*- ^ ( h > ^ ( h = V ^ ^ W ^ W 

(3.16) 

Introducing the stream function, ^ , and defining its partial derivati\'es 


as: 


9i(; 

9y 


uh 


9x 


vh 


(3.17) 


they can be substituted into the vorticity equation to give 


li . ii i_ (C±£) + it 3_ (i±£) . , ,2. + i_ . i- (i!!. 

9t 9y 9x '■ h 9x 9y '' h ^ ^ ^ 9y ^ 9x ^ph 


(3.18) 


If the Jacobian, J, is denoted as 


_ 9A 9B 

- 37 37 


Equation 3.18 can be written 


9y 9x 


(3.19) 


9c . , ,, ,C+f,^ _ . „2_ . 9 9 

— + J (i,, (—)) = A^V c + g- (^) - ^ (^) 


(3.20) 


2. Nondimensionalization 

The equations are nondiraensionalized to generalize the results, 
freeing them from the constraints of dimensionality. Nondimensional 
equations can be written in many different forms depending on the way 
in which the nondimensional parameters are formed. For example, non- 
dimensional time can be defined as: 


45 


















t 

L/U 


oo 


or 



where 


L = Characteristic length 
U = Characteristic velocity 


CO 


= Horizontal eddy diffusivity 


The choice of methods depends upon the desired result which in turn 
depends upon the way in which the equation is going to be solved. 
Physical interpretations of the parameters developed may also dictate 
the method to be used. This is because certain parameters may be 
desired which can be used to help interpret the results. 

In nondimensionalizing, reference quantities must be defined. 
The following reference quantities are the ones used in this study: 


L = Characteristic length 


= Characteristic velocity 
D = Characteristic depth 


also used in the development are: 


3 = Beta-plane approximation to the Coriolis force 


= Horizontal eddy diffusivity 


Using these quantities the following relationships between nondimensional 
and dimensional quantities are derived 


Length scale 




Depth scale 
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Velocity scale 



00 



00 


Time scale 



Nondimensional Coriolis parameter 



The remaining quantities to be nondimensionalized can be done so in terms 
of the above nondimensional quantities. They are 

Nondimensional partials of the stream function 
Knowing that 


9y 

then 


UU HD 

CO 


3Y 3y 

00 


Nondimensional vorticity 

From the definition of the vorticity 

dv _ ^ 3V ^ ^ 

^ “ 8x 3y " ^3X 3Y^ L 

thus 

2 = (^) C 


Nondimensional bottom stress 

From Appendix A the following relationship is obtained 


= lilt = 

pL h 3y HD 


(U D) 


3Y 
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therefore 


b , b 
pH ^3LU ^ ph 


Nondimensional Del operator 

The operator is defined as 
2 2 


2 2 


thus 


2 2 2 

= (L^) 

Notation is changed at this point for clarity of presentation. Sub¬ 
scripts arei now used to denote differentiation. This notation transforms 
Equation 3.20 into 

^ /5±i> . -2. 


bx bx 

CT -♦y (^) + O - V"' 

X y ^ y X 


All of the nondimensional terms are now substituted into the dimensional 
1 


equation 

U 




(U /l)z+6LF 


DH 


f + U Dfx 
X L » 


(U /L)Z+6LF 


DH 


1 

y L 


u 






3LU_ bx by 

“ l(^>y - (^>.1 


factoring yields 

U 2 U 2 

00 oo, w 

<t> - (r> ^ 


2 + 

00 

U 2 

CO 

+ (-) 

z + ei- F 

OO 

H 

H 


X 



'V = 

X 


Vo 


bx 


by 


V Z + 6U [(V“) “ ] 

" pH y pH X 


dividing by: 


U 2 

00 

<r> 


For notational purposes, small x*s and y*s will be used to denote non- 
dimensional space coordinates. 
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and defining the two nondimensional parameters 


Re 


e 


LU 

CO 

■ ' = Reynolds number 



U 

CO 


= Reciprocal of Rossby number 


the equation becomes 


- y 


^ H 


+ Y 


X 


^ = 
'' H \ 


2^ 

Re 


V^Z 


+ e 


bx 

- 


^i>y 

^pH ^x^ 


(3.21) 


Employing the definition of the Jacobian, 

7+fF 1 2 

Zt + J ^ ^ ^ ^ f Wy - 

This form is compact and neat for notational purposes but it is not 
suitable for solving numerically because it is hard to separately 
control the individual forces. Expanding the Jacobian transforms the 
equation into one which is more suitable. Expansion of the Jacobian 
yields 


J K ^ \ ^ H '^y y ^ H ^x 


H(Z4-eF) ~(Z+eF)H H(Z4-eF) -(Z4-eF)H 

= \ t -^2-- % t -^2- 


H 


H 


^ (Z+eF) (Z+eF) ‘ 

X y _ y X ^ (Z+c F) 


H 


H 


H 


. [4^ H - H ] 

2 y X X y 


(4^ Z Z ) 

'' V y X 

H ’ H ^‘x"y "y"x 


+ |- (4' F -'f F ) + [H 4' -H T 


H 


X y y X 


= I J (4',Z) + ^ J (4',F) + (^±1^) J (H,4') 


Substituting this form into Equation 3.22 























I 



(3.23) 


Zt + I J ('i'.Z) + I J (4',F) + J (H,4') = 

H 


1 2 

R? ’ " + ^ ><p->y - Wr' 


where 


= Local rate of change of vorticity with time 

J ('FjZ) = Advective term 

H 

('^>F) = Planetary vorticity tendency 

H 

J = Topographic vorticity tendency 

1 2 

V Z = Frictional term 

Re 


bx by 

^ “ (■^TT”) ] = Bottom friction 

pH y pH X 

This form is much more convenient for deleting, weighting or evaluating 
the acting forces individually. The method for incorporating the bottom 
friction will be explained in the section on the refined model and the 
complete development of the bottom friction is contained in Appendix A. 
3. Nondimensional Parameters 

Two nondimensional parameters occur in the vorticity equation. 
The first one, the Reynolds’ number, is defined as 


Re 


LU 


00 




where L, U^, and A^ have been defined previously. Physically the 
Reynolds’ number is a ratio of the viscous to the inertial forces. 
The second parameter, c, is given by 



00 
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This can be interpreted as the reciprocal of the Rossby number which is 
a measure of the effect of the Coriolis force. 

The solution of the equations depends only on the values of 
the parameters, not on the values of the individual terms in the param¬ 
eters. Nevertheless, estimates must be made of the individual terms in 
order to determine the range of values for the nondimensional parameteir 
which shou]d be investigated. Certain terms, in particular 3 and to some 

degree A^, are fairly well known. The horizontal eddy diffusivity, 

7 9 2 

has been found for turbulent flow to be between 10 and 10 cm /sec. 

7 2 

The value of 10 cm /sec will be used in this development. 3 can be 
developed from its definition. 


where 


2Q cos 9 


3 = 


^ = Angular rotation of the earth = 7.29 x 10 ^ rad/sec 


0^ = Latitude = 37 (for this study) 


r = Radius of earth = 6.38 x 10 


8 


cm 


hence, 


„ _ 2(7.29 X 10 X .8 . „„ t„-13 -1 -1 

p = -3- = 1.83 X 10 cm sec 

6.38 X 10 


Determining representative values for L and U which are the nondimensional 

oo 

length and velocity scales is not as easy. Physical measurements can be 
relied upon to give an estimate for which can be considered the free 
stream velocity or in other words the average velocity of the current 
along the coast which is unaffected by the presence of the cavity. Tlie 
best estimate that can be made indicates that a value of approximately 
one-tenth to one centimeter/second is a reasonable one. Theoretical 
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considerations are relied upon to obtain a representative value for L, 
The value chosen for L is the most questionable of all the character- 



represents the smallest possible size for any feature to be distinguish¬ 
able on the grid employed. This is because one mile is the size of one 
grid space on the arrays used to solve the equations numerically. Using 
these values, it is found that approximate values for the nondimensional 
parameters used in this study are 


-2 


Re = 1,61 X 10 


-3 


e = 4,74 X 10 


C, DEVELOPMENT OF STREAM FUNCTION EQUATION 

The vorticity equation is one equation in two unknowns; the vorticity, 
C, and the stream function, In order to solve the vorticity equation, 

another equation relating to one or both of these variables must be 
found, A second equation is derived from the definitions of the 
vorticity and stream function and equating the two, 

1, Derivation 

Respectively, the vorticity and the derivatives of the stream 
function were previously defined as 



\b = vh 

X 


^ = -uh 

y 


therefore 




X 


y 
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thus , 


hi|i - h 


C = 


XX 


X X 


hi|^ - h 


-ZZ 




finally, separating terms 




( 3-20 


Equation 3.24 can be solved for the stream function if it is converted 
to its numerical representation by substituting the proper relationships 
for the derivatives. This process is carried out in the numerical 
approach sections which occur in both the simple cavity flow problem 
and the refined model. 

2. Nondimensionalization 

Following the same procedure as in Section III-B-S and using the 
same nondiraensional terms as used for the vorticity equation, the stream 
function equation can be nondimensionalized. Substituting into the 
stream function equation yields 

z (^) - (^) y (U^LD) I (i) - ^ (V (2) + 

L HD 


4' (U D) H (f)] 
y «> y L 


collecting terms 


U 

oo 

^ <r> 


„2,„ U LD - U 

■V- ^ If H + » H 1 (-5-) 


therefore, all characteristic values cancel giving 


z = 1 (4' H + 4* H ) (3.25) 

H ^2 X X y y 

This is the desired nondimensional stream fvinction equation which is used 
in conjunction with the vorticity equation to solve for both the stream 
function and the vorticity. 
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IV. SIMPLE CAVITY FLOW 


A. GOVERNING EQUATIONS 

Tlie governing equations for the simple cavity flow problem derive 
directly from the equations developed in Section III. Two equations 
are necessciry: The vorticity transport equation. Equation 3.23, and 
the stream function equation, Equation 3.25. 


i J ('I'.Z) + ^ J (4',F) + J (H,4') = 

H 

19 bx by 


Z = 4 ^ (1' H + 4' H ) 

H X X y y 


The following assumptions are made in the simple cavity flow problem 

1. Bottom friction is negligible 

2. Depth is constant 

3. Coriolis force is negligible 

Using these assumptions, Equations 3.23 and 3.25 become 


- J (^¥,2) (4.1) 

Z = V^4' (4.2) 

where the stream function for a constant depth reduces to 
'i' = -U 

y 

^ = V 

X 

In developing the general equation, it is found that for convenience 

of form it is better to define the vorticity as 

Z = V ~ U 
X y 
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although the normal convention for the vorticity is 

Z = U - V 

y X 

Simplifying the equation eliminates the need for this unconventional 
form; therefore, the normal convention for the vorticity is used. 
Changing the convention changes the definition for the velocity in terns 
of the stream function and some of the signs in the vorticity equation. 
The partials are now defined as 

y 

= -V 

X 

Using this convention and the fact that J(A,B) = -J(B,A) the vorticity 
equation becomes 


Zt = V^z - J (Z.'i') (4.3) 

The above equation for the vorticity is in suitable form for solving 
by utilizing the Arakawa Jacobian, but a slightly different form is 
needed to solve by upwind differencing. Applying the definition of the 
Jacobian we have 


J(Z,H^) = Z - Z 'F 
X y y X 

substituting for the partials of the stream function 


J(Z,4^) = Z U + Z V 
X y 


this term can be expanded to 

UZ + VZ = (UZ) + (VZ) - Z[U + V ] 

X y X y X y 

but the bracketed term equals zero by the continuity equation, and it is 
found that the Jacobian becomes 


J(Z,'F) = (UZ)^ + (VZ)^ 


(A.A) 
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Introducing the new term into the vorticity equation, the desired form 


is obtained. 




(4.5) 


B. NUMERICAL APPROACH 

1. Upwind Differencing for Solving the Vorticity Equation 

Upv/ind differencing achieves static stability by differencing 
in the direction of the flow. Differencing of the advective terms in 
this method is always upwind of the point*in question. Considering 
Equation 4.5, forward time differencing can be used for the first term 
and the technique for representing the Laplacian can be used for the 
second. It is the advective terms in brackets which must be handled 
in a special manner. 


✓ 


- z" . 

At 


, z';,, .-zz'f. 

1 r 1 + 1,1 1 - 1,.1 1>1 

(Ax) 


+ 


...+z^^] ^-2Z^ . 
i,-|+l 1,1-1 1,1 

(Ay)^ 


- [advective terms] 


(4.6) 


where 


n is time step level 
k is iteration level in space 

The numerical scheme used to represent the advective terms depends upon 
the direction of the flow, sign of U and V, and on whether or not the 
flow reverses itself between two successive grid points. Using upwind 
differencing, the direction of the flow and occurrence of flow 
reversals must be checked at each grid point. One upwind differencing 
scheme is as follows: 
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NO FLOW REVERSALS OCCUR 


+ U 


+ V 


- U 


- V 


(UZ)^ - (UZ)^_^ 

Ix 

(VZ)^ - (VZ) ._J^ 
Ay 

(UZ),^, - (UZ). 
Ax 

(UZ).^^ - (UZ), 
Ay 


FLOW REVERSALS 


Backward differencing 
Backward differencing 


Forward differencing 
Forward differencing 


There are several methods which adequately handle flow reversals. 
The one incorporated into this study is an averaging scheme. Consider¬ 
ing the x-direction advective term, it is represented by 


Ax 


where 


^R = 


U. + U. 
1+1 1 




u. + u. , 

1 1-1 


Z„ = Z. if > 0 
R i R 


Z = Z. , if U > 0 
L 1-1 L 


\ ' "i+i \ 


Z = Z. if U < 0 
L 1 L 


The formula for the y-direction term is analogous to the x-direction 
formula. 


Interpreting the situation physically, a flow reversal appears 
as an artificial source or sink depending on the sign (- to + is a 
source of vorticity). Averaging eliminates the source or sink. 

2. Method for Solving Vorticity Equation Employing the 

Arakawa Jacobian 

This approach to the problem is very" simple. Forward time 
differencing is used for local rate of change of vorticity with time, 
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the representation for the Laplacian is used for the frictional term and 
the representation for the Jacobian is the one developed by Arakawa. 
Originally only upwind differencing was going to be used to solve the 
simple cavity flow problem but when the aspect ratio was increased to 
three the upwind differencing method would not converge. Utilizing 
the Arakawa method, solutions were obtained but at the expense of 
decreased detail. The following is the numerical equation which was 
derived for this procedure. 


- z" . 


At 


= ^ [ 


k k+l k 

. + Z. , . - 2Z. 
1+1,1 i-l»1_ij 


k k+1 k 

Z. _, + Z. . - 2Z. 


Re 


(Ax)‘ 


1 + 1,1+1 i,i-l 
(Ay)^ 


(4.7) 


10A A • (Z. . ,+Z..^ . -Z. .^,-Z.,, _,) + 

12AxAy 1+1,3 i»J-l i+l,J-l i,J+l i+l, 3 +l 

V. ... (Z.., .+Z._^T 1 --Z. 1 + 

1,3+1 1+1,3 1+1,3+1 1-1,3 1-1,3+1 

'V. T . (Z. ., 1 +Z. , ..,-Z. . ,-z. , . ,) + 

1-1,3 1,3+1 1-1,3+1 1,3-1 1-1,3-1 

'V. , (Z. .^,-Z. - .) + 

1-1,3+1 1,3+1 1-1,3 

'^i-1,3-1 ^^i-1,3 ^1,3-1^ ^ 

L Vi.j-i 

3. Optimum Overrelaxation to Solve Poisson*s Equation for the 

Stream Function 

Tlie numerical equation used to solve for the stream function is 
derived in the following manner. 

2 

Z = V ^ 

Using the representation for Laplacian 
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1 

i 

i 

' ^ 


a 




u/^ 4 - 

i"i.i 


- 2^1^ 


,k+l 

±A 


^>J+1 


2^ 


,k+l 

iJ. 


i.J 


(Ax)' 


(Ay)' 


Letting B = Ax/Ay 


(Ax)^ . 

^ J J 


Vl,j i-l,j 


2«j/k+l ^ g2 |.^k ] 

i.J i>J+l I 5 J-I 1,3 


this yields 


j,k+l 

i,j 


2(1+B^) 


['v\^ .+B^(4'^ .,,+4'^'^^ -)-(Ax)^z'? .] 

1+1,J 1~1>J 1>J+1 1>J~1 


This method for deriving an equation for the stream function is known 
as Richardson*s method. If new values are used as soon as they ‘are 
available, it becomes Leibmann’s method. A further refinement is 
achieved by using optimum overrelaxation. The formula for this 
technique is 


1,3 1,3 2(1+Bb 


- (Ax)^ Z? . - 2(1+B^)H'^ .] 
1 > J 1 > J 


(4.8) 


Ropt is the optimum overrelaxation factor which in this study has been 
found to be approximately 1.72 for the best overall results. The 
bracketed term is the residual and as the solution converges to the 
exact solution it approaches zero. 

4. Boundary Conditions 

There are many solutions to any given partial differential 
equation. Boundary and initial conditions must be specified in order 
to determine a unique solution. Figure 4.1 identifies the boundaries 
which must be specified in the simple cavity flow problem. 
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5 


1. Inflow 

2. Wall 

3. Wall 

4. Wall 

5. Lid 

6 . Outflow 


Figure 4.1 Specifications of Boundaries 


Each boundary and its associated equation are discussed separately in 
the following sections, 
a. Inflow 

(1) Velocity Profile 

To be able to specify the inflow, a velocity boundary 
layer profile must be defined. Three velocity profiles are considered 
a polynomial, logarithmic, and hyperbolic tangent. It is found that 
there is little difference in the results obtained by using these 
different profiles. The profiles are plotted for comparison in Figure 
4.2. The equations are normalized so that the values range between 
0 and the maximum velocity, 1. Over this range the polynomial and the 
logarithmic profiles agree well although the logarithmic profile does 
not approach the maximum value, 1, asymptotically as desired. 

The first profile considered is a polynomial derived 
experimentally for small scale flow in Ref. 6. It is expressed 
mathematically by the following fifth degree polynomial: 

U = .000493 + 2.408031y - 2.841989y^ + 2.796410y^ 

4 5 

- 1.904102y + .533821y 
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The second profile derives from Prandtl’s equation 
[Ref. 8]. According to Prandtl’s theory it can be shovm that the 
turbulent shearing stress becomes 

where 


K = Von Karman's constant 
u = average velocity 
Directly from Equation 4.10 



If it is assumed that the shearing stress, is constant and the 

density, p, is constant then 


therefore 


and 



constant = 


Ay. - ^ 

dy " y 


u = In y + 

To normalize this profile it is required that 


c 

II 

o 

at 

o 

II 

>> 

u = 1 

at 

y = 1 


thus 


= 1/ln (e + 1.0) 
y = ey' + 1.0 


where 


e = 2.7182818 






The velocity equation becomes 


“ ■ InleT ' O y 

The last profile considered is the hyperbolic tangent. 
This profile is considered because it is well behaved and possesses 
t\>7o desirable characteristics. First, it approaches 1 asymptotically 
and secondly near the boundary it is, to a close approximation, 
linear. In order to make the value of the profile vary from 0 to 1 
over the values of y = 0 to y = 1, the argument uy must be used. 

The resulting equation is 

U = tan h (Try) (A.] 2) 

(2) Stream Function Profile 

The stream function profile is obtained by integrating 
the velocity profile from the boundary to y. It is found that unless 
a well-smoothed stream function profile is obtained from the velocity 
profile, instabilities in the vorticity equation develop. For example, 
the following technique does not work well in deriving the stream 
function from the velocity profile. Given 



central differencing yields 


^ y 

_AJ.-1 = u 

2Ay ^ 

where in the preceding equation the subscript i = 1 denotes the first 
grid column which is the inflow boundary (Figure 4.1). Thus 


'i' 


IJ+l 


h.1-1 + 


Using this procedure, a profile like the one in Figure 4.3 is 
developed. 
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Stream Function 


Figure A.3 Unsuitable boundary layer profile 
Tills profile causes the vorticity field to alternate in sign. For this 
reason, all three velocity profiles examined are integrated analytically 
in order to obtain the stream function profiles. The three stream 
function profiles are 
POLYNOMIAL 


...... . 2.408031 2 2.8A1989 3 

^ = . 00049 3y + - - - y- - - y 


2.796410 4 1.904102 5 . .533821 6 


y - 


y + 


LOGARITHMIC 


- InVTT o ) '1" l-O) - I'O’ 


HYPERBOLIC 


T = — log (cos h Try) 


(4.13) 


(4.14) 


(4.15) 


The three profiles are plotted in Figure 4.4. It is apparent how 
similar the three profiles are. The choice of the hyperbolic profile 
for the one used in this study is based primarily on the fact that it 
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oundary layer profiles. 
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possesses the two theoretically desirable properties previously 
discussed, 

(3) Vorticity 

The development of the vorticity input equation is much 
simpler than the stream function equation. Only one assumption is 
necessary in this development, that is 


^ (^) 

3x 


= 0 


this can be approximated by 

AY.I = ill 

3x' . 3xi . 

2,j 

Remembering the relationship between the velocity and stream function 
then 


9y 




2 

1_I| 

2' 

8y i,j 


3x1 




2 

3 Y i 
2 ' 

3x l,j 


3^1 

3x^ 2,j 


Defining the vorticity as 


_ 3U . _ ^1 

‘ ‘ ^x'1.3 


substitution yields 


Z = —I 

l,j - 2 ' 


2 

3 Yi 


3y l,j 3x 2,j 
Centered differencing gives 


Z = 1 > 3+1 1 > 3-1 + 3,3 1>3 _ 2 ^ 


(Ay)' 


(Ax)' 


(4.13) 
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b. Walls 


All of the walls can be handled in a similar manner. Two 
assumptions are necessary in order to find the value of the stream 
function at the walls. 

First, it is assumed that the no slip condition applies at 

all the walls; therefore, at the walls 

U = V = 0 

'F = y = 0 
y X 

hence, ^ = constant 

Secondly, assuming that the arbitrary constant of integration is zero, 
then for all the walls 


= 0 (4.14) 

The vorticity is handled in a slightly different manner. 
Consider first the walls parallel to the X direction, at the wall it 
is obvious that V = 0, therefore 


z = ^1 

wall 9y' 


^2. 


wall 9y wall 
Expanding ^ into a Taylor’s series 


\all« -/wall ^ I 0 


this gives: 


2 ^ 


wall+1 


"wall 


(Ay)‘ 


(4.15) 


U equals zero for the walls parallel to the y direction, thus 


^wall 


^1 

9x 


= 9i!| 

2' 

wall 9x wall 


Expanding again yields: 
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wall 
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* a. 1 9 

Ax + 
wall 9x 


(Ax)' 


wall 


therefore 


^wall 


^ VaUH-l 

(Ax)^ 


(4.16) 


c- Do^TTistream Continuation 

The stream function at the outflow is determined by a simple 
extrapolation procedure. Letting IL be the outflow boundary and IL-1 
and IL-2 the two proceeding grid points, then the desired relation is 
given by 


IL-2 IL 


IL-1 


thus 


T = 2'i' - ^ 

IL IL-l,j IL-2,j 


(4.17) 


The vorticity at the outflow is assumed to be transferred 
out of the grid. This yields 


^IL ■ ^IL-1 


(4.18) 


This can be interpreted as meaning that there is no production of 
vorticity between the last two grid points, 
d. Lid 


If the Lid is considered to be a frictionless impermeable 
wall then it is implied that it is a streamline. A streamline is a line 
along which T is constant; since, the stream function is assigned a 
value at the inflow boundary then it follows that all of the other 
values along the Lid must have the same value, or 


T = Y 

i,LID 1,LID 


(4.19) 
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Determining the vorticity at the Lid is based upon two 


assumptions about the velocity at the Lid. They are 


-I 

9x' 


= 0 


LID 


9yl 


LID 


12| 

9yl 


LID-1 


from these 


^i,LID ^i,LID-1 


(4.20) 


This can be interpreted as a linear extrapolation of the U velocity 
component up to the Lid. 
e. Corners 

It is desirable to retain the effect of sharp comers in 
order to approximate reality to as close a degree as possible. The 
stream function is no problem since the value at the comer is the same 
as the rest of the wall, or zero. Tfie approach to the vorticity is a 
little more complicated and it is handled by using different values for 
the comer depending upon the point being evaluated. The scheme which 
is used is as follows 


2(¥. ... - T. .) 

A = _LiJll_ 

^ 2 
(Ay) 


Ay 




wall 


B = 


(Ax)^ 


i+l,j 


Ax 


If the vorticity at i,j+l is being evaluated the value A is used for 

Z. . when Z..., . is being evaluated the value B is used. This scheme 
I5J 1+1, J 

retains the effect of sharp comers. A similar method is employed for 
the other comer. 
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5. Convergence and Stability Criteria 


As is discussed in Section II~D, the convergence of the vorticity 
transport equation depends upon the degree of convergence of the stream 
function equation and the time step taken* The maximum allowable time 
step in turn depends upon the grid spacing, the maximum velocities and 
the Reynolds’ number. This critical time step can be determined by 
employing a discrete perturbation, stability analysis [Ref. 7]. Using 
this method it has been determined that the critical time step for the 
upwind differencing method is given by 


At 

c 


_ 1 

^ V o 

max max ^ 7 ^ 
Ax Ay Re 




(4.21) 


The time step taken must be smaller than this critical step if the 
equations are to converge. 

It is found, by making many test computer runs, that it takes 
approximately 100 relaxations of the stream function in order to allow 
the vorticity equation to converge to within acceptable limits. The 
number of relaxations is less for the method employing the Arakawa 
Jacobian than for the upwind differencing method since the former 

I 

method averages three values for the advective term, minimizing the 
growth of instabilities. 

The degree of convergence of the equations is not very strict 
since this is a preliminary study, although it is felt that the conver¬ 
gence obtained is great enough so that any further changes that might occur 
is less than the accuracy of the plotting capabilities of the Calcomp 
plotters used. In general, the equations are relaxed until the residual 
is one-tenth of one percent of its original value. The original value 
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is the one determined by the first guess field which is obtained by 
making the entire field outside of the cavity equal to the input boundary 
layer profile and all the values in the cavity equal to zero. Using 
these criteria and a time step equal to eight-tenths of the* critical 
value, it takes between twenty and forty minutes of computer time on 
the IBM 360/67 to obtain the desired results for each set of input 
parameters, which are the Reynolds’ number and the aspect ratio, 

6. Computational Sequence 

Simply stated, the procedure used to solve the numerical 
equations involves calculating the vorticity transport equation, 
relaxing the stream function equation, taking a time step and calcu¬ 
lating the new vorticity field, Tliis procedure is repeated until the 
equations are converged to within acceptable limits, A more compre¬ 
hensive list of the steps involved in this procedure is as follows: 

STEPl - Define all constants which specify the problem and/or 
minimize computational time, 

STEP2 - Specify the initial stream function by either reading in an 
initial guess field or extrapolating values from the inflow 
boundary condition, 

STEPS - (Needed only for upwind differencing method,) Specify the 

initial velocity fields utilizing the initial stream function 
field, 

STEP4 - Specify the initial vorticity field by either reading in an 

initial guess field or determining the field by using Poisson’s 
equation, 

STEPS - Calculate all the boundary conditions for the vorticity, 

STEP6 - (Needed only for upwind differencing method.) Calculate the 
velocity fields from the stream function field. 
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k+1 

STEP? “ Calculate Z using the vorticity transport equation. 

STEPS - Check the vorticity field for convergence. Exit to STEPll 
if convergence criteria is met, continue if not. 

STEP9 - Calculate the new ^ distribution using optimum overrelaxation 
holding the vorticity field constant. 

STEP10 - Advance one time step and return to STEPS. 

STEPll - Plot stream function and vorticity field computed on Calcomp 
plotters. 

7. Programming Optimization 

One of the objects of solving the simple cavity flow problem is 
to gain experience in programming problems of this nature. An important 
aspect of programming is being as efficient as possible in order to 
minimize the computer time necessary to solve the problem. Two basic 
principles to follow in order to do this are: minimize operations and 
structure the program efficiently. Minimizing operations can be divided 
into two categories. First, one can minimize computations by manipu¬ 
lating the numerical equations until they are in the most efficient form 
to be programmed. This basically breaks down to collecting like terras 
and defining new constants (to be put in the initialization part of the 
program) in order to eliminate unnecessary calculations. Secondly, all 
unnecessary computations should be eliminated. An example will demonstrate 
how this is different from the first category. The critical time step in 
the-simple cavity flow problem involves the maximum components of the 
velocity (Equation 4.21). This requires many time consuming computations 
in order to determine the maximum values on the grid. It was found that 
after five time increments the maximum values of the velocity components 
varied very little and, therefore, the time increment stabilized. It was 
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therefore unnecessary to continue these time consuming calculations 
to determine the critical time step since it remained virtually constant. 

The second basic principle of structuring the program efficiently 
can make a tremendous difference in the amount of computer time used. A 
summary of some of the things this author has found helpful is contained 
in the following list. 

1. Avoid subroutines, especially in DO loops 

2. Avoid function statements 

3. Avoid unnecessary indexing and if an index appears more 
than once in a calculation, equivalence it to a nonsubscripted variable 

4. Avoid mixed modes 

5. Avoid computed and assigned GO TO statements 

6. Convert divisions to multiplications where possible 

7. Initialize values in DATA statements 

Most of the. things in the list, which it is advisable to avoid, were 
developed for the convenience of the programmer, but the price that is 
paid for this convenience is increased computational time. 

C. RESULTS 

The complete results of the simple cavity flow problem, consisting of 
twenty-four Calcomp plots of the vorticity and stream function, is con¬ 
tained in Appendix B. These plots, twelve stream function and twelve 
vorticity, are for Aspect Ratios of 1, 2, 3, and Reynolds* numbers of 1, 
10, 100, and 1000. The plots are of only a portion of the total grid as 
sho\TO in the following diagram in which the shaded area is that area 
which is plotted. 

Contours are not plotted at equal intervals; therefore, the relative 
magnitude of the circulation within a plot cannot be determined from the 
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stream function gradients. Only the direction of the flow is indicated 
by the streamlines since the velocity vectors are tangent to the stream¬ 
lines. All plots are plotted using the same contour levels so that a 
realistic comparison between plots can be made. 



I- 4 1 -, 


PORTION OF TOTAL GRID PLOTTED 

Figure 4.5 Portion of numerical grid plotted 

The main purpose of this part of the investigation is an examination 
of the different circulation patterns which might occur. With this pur¬ 
pose in mind, four of the more interesting flow patterns are presented 
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in Figure 4.6. An indication of the effect of increasing the magnitude 
of the flow by increasing the Reynolds’ number is shown in the top two 
plots. Not only is the magnitude of the flow increased when the Reynolds’ 
number is increased, but the pattern becomes more circular and the 
center of the gyre is positioned closer to the mouth of the cavity. 'Ihe 
lower diagrams represent patterns which result if the Aspect Ratio is 
increased from 2 to 3. It is shown that the pattern starts to split in 
the case where Re = 10 and the Aspect Ratio = 3. If the Reynolds’ number 
is lowered to 1, the split becomes compJete. This latter case is inter¬ 
esting because there is some evidence that this type of pattern occurs 
in Monterey Bay. It had been speculated that the splitting of the cur¬ 
rent might be caused by the presence of the submarine canyon. But as 
seen later, this type of circulation only occurs for the simple cavity 
flow model. 

All of the vorticity plots are similar. A typical plot, for Reynolds’ 
number of 10 and Aspect Ratio of 2, is illustrated in Figure 4.7. It can 
be seen that the walls and especially the corners have a significant effect 
on the vorticity. It is doubtful that the effect of the sharp corners can 
be easily extrapolated to what occurs in the real world. 

Defining closed circulation as that which occurs in any area where the 
streamlines form a closed curve, it is interesting to note that in all 
cases closed circulation does occur. It will be seen that this differs 
markedly from the results for the refined model. The closed circulation 
varies from occurring in only the corners of the cavity to filling the 
entire cavity. If it is assumed that this simplistic model approximates 
the real world, then, important conclusions about a bay as a place within 
which to deposit effluent can be deduced. These conclusions are discussed 
in the following section. 
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Figure 4.6 Comparison of four solutions for the stream function 








































































































77 


Figure 4.7 Sample vorticity plot Re 



































D. CONCLUSIONS 


As was seen in the previous section, some form of closed circu¬ 
lation occurs in all of the solutions of the stream function obtained. 
If this is the case in reality, then, an embayment could be a poor 
place within which to deposit waste. Should the effluent be deposited 
within one of the gyres, diffusion and tidal currents would have to be 
relied upon to disperse the pollutant. Diffusion and tidal currents 
might not disperse the pollutant at a rate fast enough to keep the 
pollution below a safe level. 
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V. REFINED MODEL 


A. GOVERNING EQUATIONS 

The development of the final equations from Equations 3.23 and 3.24 
lacks only the formulation for the bottom friction. The bottom friction 
terms in the vorticity transport equation (3.23) are not in a form which 
can be numerically solved; therefore, seme other relationship must be 
found from which they can be computed. The following is a reasonable 
approach to the problem which yields a form for the bottom friction which 
can be solved numerically. 

1. Bottom Friction Equations 

Tills development, due to its length, is derived in detail in 
Appendix A. The short explanation contained here will introduce the 
method used to obtain a relationship for the bottom friction. Consider¬ 
ing Ekman’s solution for slope currents, it can be shown that 

2 

d w _ ifw 

" A A 3ri 
dz V V 

where the complex velocity field is given by 

w = u + iv 


The solution to this equation is 


ig . cosh az -v ^ 
f cosh ah dr\ 


where 


,if \ 
a - (-) 

V 


1/2 


Now, the bottom stress equals 
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T_ _ 3w I 

ph “ h 

z = 


therefore 

^ °-h" <1^' 

By integrating the velocity over depth it follows that the transport, 
Tr, equals 


Tr = (^) (|^) - (tanh ah - ah) 
r an a 

then equating 5.1 and 5.2 


(5.2) 


^ / if ') rn . tanh ah _ 

ph ^h ^tanh ah - ah 

Reducing this equation to its x^y components and nondimensionalizing 
yields 


bx 




R2 


by 


I_= F [_ R2 _ Kl ^ ] 
pH ^ H X H y^ 


R1 


(5.3) 

(5.4) 


where Rl and R2 are coefficients which are functions of ax which in 
turn is equal to 


ax 


^ , fh! 
2 


) 


1/2 


V 

Assuming that the depth at a given point, once specified, remains con¬ 
stant and f and A^ are constant, then Rl/H and R2/H need only be computed 
once at each grid point in the beginning of the numerical program and 
stored for future use in each iteration of the vorticity equation. A 
closer look at the equations shows that the bottom friction is represented 
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by depth-dependent coefficients multiplied by the components of the 
velocity. This follows directly from the fact that 

F [- R1 + R2 = - F [Rl-V + R2*U] 

n H 

The other bottom stress component takes a similar form. 

Table II contains values for R1 and R2 which are computed from 
the full equations. The frictional coefficients and approximations 
used to calculate them for large ax are plotted in Figure 5.1 and Figure 
5.2. The approximations are given in Appendix A. It can be seen from 
the plots that the approximations are excellent as long as ax is greater 
than 2.0. Below this value, the full equations must be used. Relating 
this to the depth, it is found that the approximations are good until 


h < 3/A (meters) 

V 

At very small values for the argument ax, which is dealing with depths 
on the order of centimeters, the solution oscillates about zero and 
the validity of the solution is questionable. This is not investigated, 
because at depths shallower than five meters the surf zone is encountered 
and this regime is controlled by a different set of processes. 

2. Final Form for the Vorticity and Stream Function 

Equations for the vorticity and stream function of the refined 
model were derived in Section III. These equations with slight modifi¬ 
cations, depending upon which forces are being neglected, are used in 
this section. Substituting for the bottom friction, the equations 
arrived at in Section III-l, Equations 3.23 and 3.24 become 

2^= - I J('1',Z) - I J(4',F) - 

» (5.5) 


+ eF 


[(_11 

H 


X H y'' '' H 

y 


H* - 

X 


H 
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TABLE II 


Frictional Coefficients 


ax 

R1 

R2 

.5 

-0.19995 

6.00305 

1.0 

-0.19899 

1.51133 

1.5 

-0.19509 

0.69132 

2.0 

-0.18587 

0.41530 

2.5 

-0.17090 

0.29412 

3.0 

-0.15284 

0.22951 

3.5 . 

-0.13526 

0.18872 

4.0 

-0.12014 

0.15995 

4.5 

-0.10773 

0.13849 

5.0 

-0.09756 

0.12197 

6.0 

-0.08197 

0.09836 

7.0 

-0.07059 

0.08235 

8.0 

-0.06195 

0.07080 

9.0 

-0.05517 

0.06207 

10.0 

-0.04972 

0.05525 

15.0 

-0.03325 

0.03563 

20.0 

-0.02497 

0.02628 

25.0 

-0.01998 

0.02082 

30.0 

-0.01666 

0.01723 

35.0 

-0.01428 

0.01470 

40.0 

-0.01250 

0.01282 

45.0 

-0.01111 

0.01136 

50.0 

-0.01000 

0.01020 

100.0 

-0.00500 

0.00505 

200.0 

-0.00250 

0.00251 

300.0 

-0.00167 

0.00167 

400.0 

-0.00125 

0.00125 

500.0 

-0.00100 

0.00100 

1000.0 

-0.00050 

0.00050 

2000.0 

-0.00025 

0.00025 
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Figure 5.1 Frictional coefficient, R1 
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Figure 5.2 Frictional coefficient, R2 
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and 


Z 


^ V^Z - ^ ('1' H -f- f H ) 
H „2 X X y y 


(5.6) 


Taking the indicated partials in the bottom friction term gives 


Zt = — V^Z - I J(f ,Z) - I JCH*,?) - 

H 


(5.7) 


+ eF [|^ + J(|i,H') +(1^)4' + (|^) 4'^] 

y ^ X 


The equations for the refined model contain the volume transport 
stream function rather than the velocity stream function which is con¬ 
sidered in the simple cavity flow problem. This study is interested 
only in the direction of the flow and net the volume transported. The 
velocity stream function has the property that the direction of the flow 
is at all points tangent to the streamlines. Here^ it will be proven 
that the volume transport stream function has the same property. 
Considering 


V-V¥ 


uf+v|I 

dx oy 


therefore 


V-VH^ = U (+ UH) + V(-IJH) 


thus 


V-V'F = 0 


This means that the velocity is perpendicular to the gradient of the 
volume transport stream function which in turn means the velocity is 
everywhere tangent to the streamlines, 

3. Convergence and Stability 

The rate of convergence for the refined model, as expected, is 
much slower than for the simple cavity flow problem. To obtain 
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convergence to within one percent of the original residual, it takes 
between two and three hours of computer time. On the average, the 
equations are converged until the residual is one-half of one percent 
of the original value. The proper time £;tep is determined experimentally 
and is found to be about one-fiftieth of the time step used in the 
simple cavity flow section. Convergence to within one percent of the 
original residual required approximately 600 time steps. 

Although the Arakawa Jacobian is a stable representation, it is 
found that the bottom topography introduced must be smoothed, eliminating 
some detail, in order to obtain the desired degree of convergence. 


B. INVESTIGATION OF THE RELATIVE IMPORT.^NCE OF THE CORIOLIS FORCE AND 
BOTTOM FRICTION 


One question which occurred during the development of this study v^as: 
how importcint are the Coriolis force and the bottom friction? To inv6;s- 
tigate this question the nondimensional parameter, e, and the non- 
dimensional. Coriolis force, F, is examined since both forces contain 
these termf;. 

1. Coriolis Force 

The magnitude of the Coriolis force and, to a degree, the bottom 
friction depend upon the magnitude of the nondimensional parameter, e. 
Recalling that 


e 



(5.8) 


where 


and 


3 = 


20 cos 9 

o 

r 


F = F + Y 
o 
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where 


Y -= Direction of true north 

F ^ tan 9 
o L o 

0, " (7.29x10 ^ rad/sec) angular rotation of earth 

g 

r == (6.38x10 cm) radius of earth 

L == (1.609x10^ cm) characteristic length = one mesh length 
in X direction 

9 = (37^) latitude for this study 

o 

Carrying out the indicated arithmetic operations yields 

-13 -1 -1 

8 = 1.8 X 10 cm sec 

and 


which reduces to 


F = 


F = 


6.38 X 10^ 

1.61 X 10^ 

3 X 10^ + Y 


tan 37° 


+ Y 


It can be seen from 5.8 that once the characteristic length is deter¬ 
mined, then e depends only upon the characteristic velocity, U^. Fron 
the definition for the Reynolds’ number, a relationship for can be 
determined. 


ReA^ 


Substituting into 5.8 gives 



Assuming a value of 10 for and substituting the values obtained for 
8 and L, then e becomes 


G 


7.64 X 10 
Re 
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This is a convenient relationship because solutions for various 
Reynolds’ numbers are desired; therefore, once the Reynolds’ number 
is specified, the value for c is determined. 

Assuming, at first, that the Coriolis force is constant, the 
term appeairing as a coefficient in front of both the bottom friction 
and topographic tendency terms is 

cF = cF 

o 

Direct substitution yields 


The inclusion in the model of the constant term, cF, is necessary and 
presents no problems. Now, the effect of the inclusion of the variati.on 
of the Coriolis term with Latitude is examined. The term being examined 
is given by 


f f [4' F - 'V T ] 

H H X y y X 

It has been determined that 


eF = 3 X 10^£ + Ye 

Thus, if it is assumed that the grid is aligned such that Y is in the 
direction of True North, then 


(eF) 


y 


7.64 X 10 

Re 


but 


(eF) = eF 

y y 


therefore 


■| [T F 
H X y 



7.64 X 10 lx 
Re H 
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This term is at least several orders of magnitude less than the other 
terms in the equation and is, therefore, negligible. Thus, although 
the variation of the Coriolis force is m^gligible due to the scale used 
in this study, the presence of a constanc Coriolis force will be shown 
to have an effect upon the solutions which are obtained. 

2. Bottom Friction 

Examining the bottom friction terms in Equation 5.7 it can be 
seen that the frictional coefficients, RL and R2, appear in all of thei 
expressions. Looking at Table II it is seen that the frictional coeffi¬ 
cients rapidly become small and the change of the coefficients iii the 
horizontal direction, due to the change in the bottom topography, becomes 
negligible. Considering cF, with Reynolds’ numbers greater than one, 
this factor which is multiplied into the frictional term is less than 
.2. Thus at Reynolds’ numbers greater than one or at depths greater 
than thirty meters the bottom friction becomes negligible. Since this; 
model uses a grid spacing of 0.5 miles in the Y direction and 1.0 miles 
in the X direction, the area covered and the topographic changes are 
large. The only place where the bottom friction would have any signifi¬ 
cance is at the few grid points immediately surrounding the boundaries 
where conditions are only approximated. At very shallow depths, less 
than five meters, the equations no longer represent the current regime 
because the surf zone is encountered. Also, considering that the 
velocity profile, irrespective of bottom friction, is generally a maximum 
at the surface and decreases with depth, the effect of bottom friction 
is minimized in this situation compared to the situation where the 
velocity is constant with depth. For these reasons, it is concluded 
that it is unrealistic to assume that the inclusion of the bottom 
friction would increase the accuracy of this model. 
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C. NUMERICAL APPROACH 


1, Method for Solving the Vorticity Equation 

The full vorticity equation, Equation 5.7, neglecting the 
variation of the Coriolis force with Latitude and the bottom friction, 
reduces to 


V^z - j(h,H') 

H 


(5.9) 


in finite difference form this becomes 

.„k+l „k 

+ 


k+1 k 

„n+l _n . Z7._ .+Z^ I .-Tir, . 

AT Re ^ ,, .2 

(Ax) 


• ,-2Z^ , 

i>3+l i>J-l i,J 


(Ay)- 


(5.10) 


, (Z. . + eF) 

- J(T,Z) + -J(H,'i')] 

i 5 j H. . 

where the Jacobians in the above equation are solved by Arak^wa’s 
technique v^hich is outlined in Section II-B-3. For programming purpo£;es, 
the above equation is manipulated to yield a similar equation which ml.ni- 
mizes number of computations needed to obtain a solution. 

2. Optimum Overrelaxation for Solving the Stream Function Equatio n 

Substituting the appropriate finite difference forms into 
Equation 5.6 yields 


Z = 


H. 


[■ 


1 + 1,3 i - l >1 ■ ■ --11 


i + i.-i+l 1,3-1 i,.i 


i>3 


(Ax)- 


] 


(Ay)- 


1 (T.-t . .)(H.,, .-H. , .) ('V. . ,)(H. . ,) 

A _ r 1+1,3 1-1>3 1+1,3 1-1.3 . 1.3+1 1»3-1 1,3+1 _i 

n L n ^ o J 

H. . 4 (Ax) 4 (Ay) 

1 jj 


An expression for .is desired; therefore, first multiply by (Ax)' 

^ >3 

and then collect terms obtaining: 
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2 2 

2 [1 + (7^) ] .) + ( t ^) (H*. . . 1) 

Ay 1,J 1+1, J 1-1 ,J Ay 1,J+1 1,J-1 

- (Ax)^ H, .Z? . - I 

"l,j 4 (Ax>2 


a. ..,-'1'. . ,) (H. , ,) 

+ . ■?- > .3.. +l ...i^J--l_ 1,3+1 1,3-1 

4(Ay)^ 


] 


AX 

Letting B = — and employing Optimum Ovei*relaxation 


^k+l 

i,j 


, V k Root k k+1 2 k k+1 

(1.0 - Ropt) 4^. . Pi/ + B^(>i/ + '!/. ) 

" " ,2^ ^ i+lsj 1"1>J i>J+l 


2(1+B ) 


- (Ax)^ H. .z” . 

i,J i,J 


. ^ k+1 ,, 

.-'y. 1 .-H. ^ .) 

1+1,3 1-1,1 1+1,3 1-1,3 

4H, , 

1,1 


k k+1 ^ , 

.^1 - 'I'. . 1 ) (H. 


- H• . 1 )] 


4H. , i,3+l i,j-l i,j+l 1,3-1' 

) J 


Defining the constants and reorganizing in order to minimize computations, 
this equation is used to calculate the stream function. 

3. Boundary Conditions 
a. Inflow 

In this model the stream function is not related to the 
velocity, U, as in the simple cavity flow but to the volume transport, 

UH; therefore, the stream function is a function of both the depth and 
the velocity. It is assumed that the inflow boundary layer is completely 
dominated by the bottom profile with frictional effects being negligible. 
This being the case, the input velocity is assumed constant and the stream 
function boi.indary layer is determined directly from the bottom topography. 
Knowing that 

^ = --UH 

y 
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then one procedure for determining the input boundary layer is 



and 



This procedure for determining the input boundary layer is found to be 
unsatisfactory as it is in the case of the simple cavity flow. To 
solve the problem of specifying the inflow stream function boundary 
layer, the grid is extended five grid spaces upstream and four spaces 
downstream as transitional regions. The bottom is specified as linear 
at the artificial inflow boundary gradually changing into the "true” 
bottom profile. When this is done, a well behaved input profile can be 
determined by 


T = -UH 

y 


thus 



UH dy 


y, 


o 


but U is constant and H = CY, therefore 



Y dy 


y, 


o 


finally 



2 

f = -uc 



where C is the slope of the artificial inflow bottom 
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It follows directly from Equation 3.11 that the vorticity 
is zero at the inflow since the V component of the velocity is zero and 
the U component is constant. 

b. Lateral Boundaries 

The rectangular shape of the cavity is retained in the 
refined model. It is assumed that the effect of the bottom topography 
will simulate the effect of realistic boundaries. This is achieved by 
using the actual coded bottom but specifying a minimum depth to be used 
whenever the depth becomes too shallow or the boundary goes overland. 
This minimum depth is used until the wall is reached. Diagramatically, 
this is illustrated in Figure 5.3. Using this method, the problems of 
programming variable boundaries and its associated large increase in 
computer time is bypassed. In doing this the boundary condition at the 
wall still must be specified. 


The stream function at the v^alls is still equal to zero since 


no slip conditions are applied wiiich implies that the velocity at the 
walls is zero. The vorticity can be determined in the following manner. 
First, consider walls alligned in the X direction. At the wall, V = 0 
and 



thus 



Z 


'wall 


y 


H 


and 


Z 11 = ^ 

wall H yy 


Considering Taylor’s expansion: 
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Figure 5.3 Grid for refined model. 
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4' 

wall+1 




Ay + -r- 4' 

2 yy 


I (Ay)^ 

'wall 


therefore 

>''' (4y)^ 


Substitution yields 

z 

wall H. . 


^^wall-Hl 

(Ay)^ 


Vertical walls are handled in a similar manner yielding 


\'all 


2^ 

wall+1 
1,3 (Ax) 


Sharp corners are handled in the same way they are handled 
in the sim})le cavity flow problem. It is expected that the effects of 
the boundaries on the solution will be diminished due to the effect of 
the bottom topography, 

c. Downstream Continuation 

The same method as the one employed for the simple cavity 
flow problem is used. This method gives 

w = w = 2 ^ - ¥ 

outflow IL,J IL-1,J IL-2,J 

and 

^outflow " ^IL,J " ^IL-1,J 

d. "Lid" 

It is assumed that the *’Lid” is far enough from the cavity 
so that the assumption that there is no volume transport across this 
imaginary boundary will not seriously affect the flow patterns created. 

It is assumed that the current at the ”Lid’^ is constant and parallel to the 
”Lid'*. This makes the **Lid'* a stream line and as in the simple cavity flow 
problem: 
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^i,LID ^i,LID-1 

e. Bottom Topography 

Two basic bottom topographies are used in this study. The 
first bottom topography (Figure 5.4) used is a hypothetical one which 
is very simple in order to minimize any instability it might generate. 
This bottom is obtained by making the depth at the ”Lid" a constant 1600 
meters and assuming a linear decrease of depth along the inflow boundary 
to zero at the first grid after the wall. This makes the value at the 
wall the minimum non-zero depth. All of the other columns exterior to 
the cavity are equivalenced to the inflow column. The entire grid is 
then relaxed, holding the boundaries constant until a steady state is 
reached. The presence of the cavity created the bottom topography in 
Figure 5.4. 


To obtain a bottom topography which simulates that of 
Monterey Bay, a bathymetric chart of the bay is overlain by a 50 x 80 
grid and the depths are coded at each grid intersection and then placed 
on computer cards. The field is then read into the computer and stored 
for use in the program as needed. The maximum depth of 1600 meters is 
chosen because it is approximately the maximum level of no motion used 
in geostrophic current calculations. The field is smoothed by averaging 
the surrounding four grid points in order to obtain the central grid 
value, i.e.: 


H = 1+1 > 3_IzIiO_ 1,3+1 1,3-1 

i,3 ^ 

It is necessary to smooth the field in order to eliminate inconsistencies 
in the bottom topography. In order to retain the effect of the canyon, 
the boundaries of the 1600 meter contour are preserved. This boundary 
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Figure 5.4 Hypothetical bottom topography 

























































is marked by the dark line in Figure 5,5. Figures 5.6—5.9 show a 
progression from the original field to one smoothed five times. It 


can be seen that after only five smoothing passes only the general 


outline of the canyon is still preserved. 


4. Computational Sequence 


The computational sequence for this portion of the study is v(iry 


similar to the one used in the simple cavity flow section with a few 
additions. The first difference is that the bottom topography must 
either be generated within the program or must be read into the program 
from some external source: tape, disc, cards, etc. The other main 
difference is that the vorticity transport equation is much more compl'.ex 
and some of the dimensioned variables occur more than once in the 
equation. In order to minimize computer time, at each time step, thesje 
multiple-subscripted variables are equivalenced to non-subscripted 
variables before the vorticity transport equation is solved. 

D. RESULTS 

The vorticity equation used for the refined model is given by 





(5.11) 


H 


In the above equation the stream function is the volume transport stream 
function. There are two cases which develop naturally from this equation. 
Case I occurs when cF equals zero. Physically, this is the situation 
that arises when the Coriolis force is neglected. When cF is assumed to 
be a constant other than zero, then only the change of the Coriolis force 
with latitude is being neglected. This is case II. Solutions for both 
cases are obtained and discussed in the next two sections. A measure of 
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Figure 5.5 Bottom topography simulating Monterey Bay. 


























































100 


Figure 5.6 Bottom topography with no smoothing passes. 
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Figure 5.7 Bottom topography with one smoothing pass. 
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Figure 5.8 Bottom topography with two smoothing passes. 
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Figure 5.9 Bottom topography with five smoothing passes. 
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the importance of the various terms in Ecuation 5.11 can be obtained by 
varying the Reynolds’ number and cF. 

1. Neglecting the Coriolis Force 

For this case, results are obtained for both bottom topographies 
discussed in Section V-'S-e. In the solutions obtained here the flow is 
from the left to the right representing flow from the South to the North. 
For the hypothetical bottom topography illustrated in Figure 5.3, due to 
its symmetry, the solutions are the same for flow in either direction. 

This is true only because the input velocity is assumed constant and 
therefore the input vorticity is zero. This is tested in several computer 
runs by reversing the position of the inflow and outflow and differencing 
in the opposite direction. No differenctiis in the solutions develop from 
reversing the direction of the flow, but it is expected that if gyres 
occur v;hich are not in the center of the cavity then the two solutions 
would have the gyre positioned on opposite sides of the cavity. If there 
is an input velocity boundary layer, different solutions would be obtained 
for flows in opposite directions because of the change in the sign of the 
input vorticity. 

Using the hypothetical bottom topography illustrated in Figure 
5.3, the volume transport stream function plot in Figure 5.10 is the 
result obtained. Computer runs for Reynolds’ numbers of .01, .1, 1, 10, 

100 are made and all of the results are qualitatively identical to this 
plot. The plot of the volume transport stream function obtained when 
the bottom topography simulating Monterey Bay is used is contained in 
Figure 5.11. The results for this bottom topography are also identical 
over the range of Reynolds’ numbers from .01 to 100. 

It can be seen from the results for the hypothetical bottom 
topography that within the bay, the streamlines follow the bottom contours. 
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Figure 5.10 Refined model, hypotlietical lioLtoin topograpliy, volume transport stream function 
Re = .01-100 eF = 0. 
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Figure 5.11 Refined model, topography simulating Monterey Rey, volume tianspori. Suj-eam j.unci.xon 
Re = .01-100 £F = 0. 
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This means that the direction of the flow within the bay is generally 
along the bottom contours; but at the larger depths encountered outside 
of the bay, the direction of the flow does not follow the contours. The 
fact that the flow does not follow the contours outside of the bay while 
it does inside the bay is undoubtedly, not only, the result of shallower 
depths within the bay, but also, due to the fact that the volume trans¬ 
port in the bay is much slower than outside of the bay. For the bottom 
topography simulating Monterey Bay the entire cavity is filled with a 
gyre. Again, no significant difference is discerned between the solutions 
obtained o\er the range of values investigated. It is also noticed that 
the cavity seems to cause a disturbance downstream from the cavity. 
Although the field dov/nstream is not long enough to be absolutely certain, 
it is reasonable to assume that this is a wave induced into the 
flow by the bottom topography. For comparison purposes, the simple 
cavity flow problem is solved for the case of Reynolds* number equal to 
.01 and As])ect Ratio of 2. The result is contained in Figure 5.12. It 
is evident what a significant difference the inclusion of the bottom 
topography has upon the resultant solution. 

2. Assuming a Constant Coriolis Force 

After examining the case where the Coriolis force is neglected, 
the effect of including a constant Coriolis force is investigated. A 
range of values of the Reynolds’ number and Coriolis force must be 
investigated because these parameters depend upon the characteristic 
scales used. Since some of the scales, such as the length scale, cannot 
be determined exactly, solutions must be obtained for a range of values 
of the nondimensional parameters. It is then assumed that the solution, 
which models the true circulation the best, falls within this range of 
values. 
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Figure 5.12 Simple cavity flow, stream function, Re 




































































Using the hypothetical bottom topography shown in Figure 5.3, 
solutions are obtained for eF equal to .23 and 23. These are contained 
in Figures 5.13 and 5.14, respectively. Each figure represents solutions 
for Reynolds* numbers of .01 and 1 since the solutions do not differ over 
this range of values. As can be seen a change of two orders of magnitude 
of cF has some effect on the flow pattern generated; whereas, a change 
of two orders of magnitude of the Reynolds* number has a negligible 
effect upon the solution. This means that over this range of values the 
effect of the Coriolis force interacting with the bottom topography, 

H 

has a larger effect upon the solution than the frictional term. 



Still, no closed circulation occurs in any of the solutions obtained for 
the hypoth(itical bottom topography. 

Th(i next step is to obtain solutions for various values of the 
Reynolds* number and the Coriolis force for the bottom topography simu¬ 
lating Monterey Bay, Figure 5.4. Solutions are obtained for the follow¬ 
ing values of the parameters: 


eF 

M 

Figure 

-2 



.23 X 10 

.01, 1 

5.15 

.23 X 10^ 

1—1 

I—1 

O 

5.16 

1.00 X 10^ 

1—1 

1—1 
o 

5.17 


It is immediately obvious, that with the Monterey Bay bottom topography, 
closed circulation occurs for all of the solutions obtained. Also, the 
solution for the hypothetical bottom fails to yield closed circulation 
over the same range of values for the two nondimensional parameters, Re 
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Figure 5.13 Refined model, hypothetical bottom topography, volume transport stream fvinction 
Re = .01-1.00 eF = .23. 
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Figure 5.14 Refined model, hypothetical bottom topography, volume transport stream function 
Re = .01-1.0 eF = 23.0. 
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function. Re = .01 cF = .23 x 10~^. 
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Figure 5.16 Refined model, bottom topography^simulatlng Monterey Bay, volume transport stream 
function. Re = .01 eF = .23 x 10^. 
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Figure 5.17 Refined model, bottom topography simulating Monterey Bay, volume transport stream 
function. Re = .01-1.0 eF = 1.0 x 10^. 



























































and eF. It can be concluded, then, that the closed circulation is the 
result of the presence of the Monterey Submarine Canyon and that the 
bottom topography has the most significant effect upon the circulation 
pattern generated. 

E. CONCLISIONS 

The foundation for the development of a good model for investigating 
the current-driven circulation in an embayment has been established. 
Included in this model are the effects of advection, planetary vorticity 
tendency, topographic vorticity tendency and lateral shear stress. The 
relative order of magnitude of the bottom friction and the change of 
Coriolis force with latitude is investigated and determined to be 
negligible! compared to the other acting forces and assumptions that are 
made in the development of this first stage of the model. 

There are two areas in particular where the model could be improved. 
First, effects of the inclusion of the bottom friction should be investi¬ 
gated further. Effects of the bottom friction diminishes rapidly with 
increasing depth, but in shallow areas around the boundaries and within 
the bay, the bottom friction might have a significant influence upon the 
circulation patterns. Second, some of the boundary conditions should be 
refined. The walls representing the coastline might be represented in 
a more realistic manner. This would go hand in hand with the inclusion 
of the bottom friction. Also, the inflow and outflow should be represented 
by less restrictive equations. 

The representation for the bottom friction is an approach suggested by 
Dr. Jerry Galt and developed in this thesis by the author. It is believed 
to be a new approach and one worth investigating further. This approach 
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is of special interest because the bottom friction is a linear function of 
the velocity making it amenable to analytic solutions. 

It is assumed in this model that the input vorticity is zero; because 
of this, the solution does not depend upon the sign of the vorticity gen¬ 
erated at the boundary which is a function of the direction of the flow. 

If a velocity boundary layer exists, as many studies propose, then the 
vorticity will not be zero at the inflow and the solution will, in addition 
to the other variables, become a function of the direction of the flow. 

This is a case which should be investigated not only to find out the effect 
of a boundary layer upon the solution but to determine the effect on the 
solution of changing the direction of the flow. 

Pvesults of this study indicate that the occurrence of closed circu¬ 
lation in !'onterey Bay is a distinct possibility. This investigation 
indicates that the presence of the submarine canyon seems to increase the 
possibility of the occurrence of closed circulation. 

In conclusion, it can be stated that this study indicates that the 
bottom topography is the controlling factor in determining what type of 
circulation pattern is to be expected in a given bay. Over the ranges of 
values examined for Re and cF, the change of the vorticity from one case 
to the next had little effect upon the solution of the stream function. 

This is because in deeper v/ater the volume transport becomes much larger 
than the vorticity and, therefore, the stream function equation becomes 

i (Y H + 'F H ) = 0 

H ^ ^ y y 

which is, for all practical purposes, independent of the vorticity., This 
puts additional constraints upon the types of flow patterns which may occur. 
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APPENDIX A 


DEVELOPMENT OF BOTTOM FRICTION 


Assuming that the coefficient of vertical eddy diffusivity, A^, is 
constant, the bottom stress is proportional to the change of the mean 
velocity gradient with depth. Using Ekman*s solution for slope currents 
an expression for the velocity as a function of depth can be obtained, but 
in this relationship the velocity is also a function of the slope of the 
sea surface. In order to eliminate the sea surface slope, which is not 
easily obtainable, the velocity can be integrated over depth to give an 
expression for the volume transport. Equating these two relationships 
an expression for the bottom stress which does not depend upon the slope 
of the sea surface can be obtained. The first step in this development 
is to obtain the analytic solution for the velocity from the following 
three equations 


fv + A 

V 



1 

p '6x 


-fu + A 

V 



i ^ 

P 


(A-1) 


(A-2) 


P = “pgh (hydrostatic pressure) (A-3) 

These equations can be derived from equations II~1 and 11-2 with the 
addition of the following assumptions 

1. Steady state conditions prevail 

2. Non-linear terms are negligible 

3. Lateral turbulent stress is negligible 

Ekman solved these equations analytically in the following manner: 
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and 


3P 3h 

3x 3x 


3P 3h 

3y 3y 


Substitution yields 


i!u 

3 z 2 


= _ ^ 
A A 9x 

V V 


thus 


defining 


it follows 


3 V _ Ju 

. .2 ^ A A 3y 

3z V V 


(u+iv) = (u+iv) + -f— 4- i -^) 


3z 


V 


w = u + iv 


3^ 3h ^ . 3h 


A 3x 3y' 

V 


3y 


d w ^ fw g 8C 

,^2 A A 3n 

dz V V 


Ekman found that the solution to this equation is 


(A-4) 


„ = _ /M-) f Cosh az _ H 
cosh ah 3n 


where 


a = (^) 

V 


1/2 


(A-5) 


Substituting A-5 into equation A-4 readily shows that it is a solution. 

Given an expression for the velocity, expressions for the bottom 
stress and the transport can now be found. The bottom stress is given by 
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Ph-h 


since 


9w 


therefore 


3Z = ^ ^ 


^ = tanhah|^ 


(A-6) 


Determining the transport is as easy as the determination of the bottom 

friction. Knowing that 

o 


then 


thus 


= J wd: 


-h 


■■'r ■ / > < 


cosh aZ 
cosh ah 


- 1) dz 


-h 


T„ = (¥) (If) [ 


r f 9n cosh ah a 


1 1 • u 

sinh 32 


- Z 


-h 


-h 


finally 

Equating A-6 and A-7 gives 


= C-^) (^) (—) (tanh ah - ah) 


(A-7) 


r \' ir 11 - 1^1 


from this the desired expression for the bottom friction is obtained 


^ A , 

^ V 2 ^ r 

‘T' ^ ^ T [ 

ph K - ^ 


tanh ah , 
r ‘'tanh ah -- ah^ 


(A-8) 


This form is not convenient for programming. Transforming equation A-8 
into a more convenient form involves separating it into its x,y components, 
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This is accomplished in the following manner: first, the transport is 
defined in terms of the stream function 


T = T H- i T 
r X y 


but 


T = m = 

X y 


T = VH = Y 

y X 


thus 


secondly, the 
Defining 


T = + i 'F 

^ y X 


hyperbolic terras must be separated into its components, 


9 1/2 

, .1/2 .fh^. .1/2 

ah = 1 —) =1 X 

V 


therefore 


b 


T 




1/2 

TANK i ' X 


. 1/2 . 1/2 ^ 

TANH 1 X - 1 X 


(A-9) 


where 


il/2 


COS 45 


o 


+ i sin 45 


o 


‘^4. • 4.- 

= — + ^~2 


Substituting for the arguments of the hyperbolic terms yields 

1/2 

TANH i ^ X ^ TANH (ax+iax) _ 

.1/2 .1/2 TANH (ax+iax) - (ax+iax) 

TANH 1 X - 1 X 

through fundamental hyperbolic trigonometric identities and some simple 
complex variable manipulations the following relationship can be arrived 
at 


TANH (ax+iax) _ = R1 + i R2 

TANH(ax +iax) - (ax+iax) 


in which: 
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R1 = 


R2 = 


G-C 


G + 2 [A-C] 
D 

G + 2 [A'-Cf 


and 


A = ax [COSH 2ax + COS 2ax] 
C = SINH 2ax + SIN 2ax 
D = SINH 2ax SIN 2ax 


G = 


(COSH 2ax - COS 2ax) 

ax ■ 


substituting into equation A-9 gives 


(A-10) 


^ (_vj» +14/ ) (Ri + 1 R2) 
ph h y X 

expanding terms yields 

bx by ^ 

+ i ^ [-R1H' + R24' + i (-R2'F - R14' ) ] 

ph ph h X y X y 

equating real and imaginary terms leads to the desired bottom frict 

components, i•e.: 


ion 


bx 

V- = ^ (-Rl'{' + R2'f' ) 

ph h X y 


and 


f 

-V- = r (-R2'}' - R14' ) 

ph h X y 

Approximations for the frictional coefficients, Rl and R2, can be derived 
when 2ax is greater than 4 by eliminating negligible terms in the 
equations A-10, The approximations can be shown to be 

1 ~ ax 


Rl - 


2 (ax) 2ax + 1 
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ax 


2 

2 (ax) - 2ax + 1 

After these equations were developed it was determined from the 
relative order of magnitude of the terms involved that the bottom 
friction is negligible; therefore, this representation was not tested. 

Plots of the frictional coefficients and tables of it’s values are 
contained in Section V. By way of comparison, though, it can be shown 
that the results from this method can be interpreted to agree with the values 
for the bottom friction used by Hansen^ who represented the bottom friction 
by 

b / / 2 ^ 2 

T = (ru vu +v , rv vu +v ) 


where 

r = 3 X 10 

This representation is the same as 


-3 


T = 3 X 10 X V 

where V is the velocity vector. The method described in this thesis gives 
approximately the same value for *’r" over the depth range of 2000 - 3000 
meters if it is assumed that A = 100. A big difference is that the 

V 

2 

coefficient in this study is multiplied by V instead of V , but if it is 
assumed that the mean currents over depth investigated are of the order of 
magnitude of 1 cm/sec then the results from both methods are similar. It 
is impossible to compare the two exactly because the methods for computing 
the bottom friction are dissimilar, Tliis method does allow the bottom 
friction to be represented in a rational manner and has the added 
attraction of being linear in velocity making it amenable to inclusion 
in analytical solutions. 


Hjalter Hansen, Hydrodynamical Methods Applied to Oceanographic Problems , 
Proceedings of the Symposium on Mathematical-Hydrodynamical Methods of 
Physical Oceanography, p. 25-34. Institut fur Meereskunde der Universitat 
Hamburg (1962). 
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APPENDIX B 


SOLUTIONS FOR THE SI>IPLE CAVITY FLOW PROBLEM 
This appendix contains Calcomp plots of solutions to the simple cavity 
flow problem. There are twenty--four plots, twelve of the stream function 
and twelve of the vorticity. Two different techniques are used to obtain 
these solutions. Upwind differencing is used to obtain as many solutions 
as possible but when this technique failed to converge the second method, 
employing the Aralcawa Jacobian for solving the advective terms, is used. 
The following is a breakdown of the technique used to obtain each 


solution. 



Aspect Ratio 

Reynolds’ Number 

Technique 

1 

1 ,10,100,1000 

Upwind differencing 

2 

1 ,10,100 

Upwind differencing 

2 

1000 

Arakawa 


3 


1,10,100,1000 


Arakawa 
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Figure B.l Simple cavity flow, stream function, Re 
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Figure B.2 Simple cavity flow, vorticity, Re 
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Figure B.3 Simple cavity flow, stream function. Re 
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Figure B.4 Simple cavity flow, vorticity, Re 
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Figure B.5 Simple cavity flow, stream function, Re 
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Figure B.6 Simple cavity flow, vorticity, Re 
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Figure B.7 Simple cavity flow, stream function. Re 
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Figure B.8 Simple cavity flow, vorticity, Re 
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Figure B.9 Simple cavity flow, stream function, Re 
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Figure B.IO Simple cavity flow, vorticity, Re 
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Figure B.ll Simple cavity flow, stream function, Re 
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Figure B.12 Simple cavity flow, vorticity. Re 
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Figure B.13 Simple cavity flow, stream function. Re = 100, 
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Figure 3.14 Simple cavity flow, vorticity, Re = 100, 
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Figure B.15 Simple cavity flow, stream function, Re = 100, 
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Figure B.16 Simple cavity flow, vorticity. Re = 100, 
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Figure B.17 Simple cavity flow, stream function, Re = 100 
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Figure B.18 Simple cavity flow, vorticity. Re = 100, 
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Figure B.19 Simple cavity flow, stream function. Re = 1000, 
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Figure B.20 Simple cavity flow, vorticity, Re = 1000, 
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Figure B.21 Simple cavity flow, stream function, Re = 1000 
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Figure B.22 Simple cavity flow, vorticity. Re = 1000, 












































CO 

il 

pq 


146 


Figure B.23 Simple cavity flow, stream function. Re = 1000, 
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Figure B.24 Simple cavity flow, vorticity, Re = 1000, 
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